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ABSTRACT 


This  thesis  reviews  two  algorithms,  one  for  the 
computation  of  a  Reflexive  Generalized  Inverse  and  the 
other  for  the  computation  of  the  Pseudoinverse  of  an 
arbitrary  matrix.  Upper  bounds  are  obtained  for  the 
round-off  error  incurred  in  performing  these  computations. 
Several  Reflexive  Generalized  Inverses  and  Pseudoinverses 
are  computed,  and  a  comparison  is  made  between  the  pre¬ 
dicted  bounds  and  the  actual  round-off  errors. 


. 


ACKNOWLEDGEMENTS 


I  express  my  appreciation  to  Professor  S.  Charmonman 
for  the  guidance  given  me  in  the  preparation  of  this  thesis, 
to  Professor  U.M.  von  Maydell  for  her  interest  and  assistance 
in  this  topic,  and  to  Professor  D.B.  Scott,  Head  of  the 
Department  of  Computing  Science,  for  providing  computing 
facilities  and  financial  assistance  while  this  research  was 
being  done.  I  also  wish  to  thank  the  National  Research 
Council  of  Canada  for  financial  assistance  for  carrying  out 
this  research. 


TABLE  OF  CONTENTS 


Page 

CHAPTER  I  -  INTRODUCTION 


1.1  History  1 

1.2  Purpose  of  the  Study  4 

CHAPTER  II  -  BASIC  THEORY 

2.1  Matrix  Algebra  5 

2.1.1  Notation  and  Basic  Definitions  5 

2.1.2  Permutation  Matrices  6 

2.1.3  Hermite  Normal  Form  7 

2.1.4  Matrix  Norms  9 

2.2  The  Elimination  Methods  of  Gauss  and  Jordan  10 

2.2.1  Gaussian  Elimination  10 

2. 2. 1.1  Partial  Pivoting  12 

2. 2. 1.2  Complete  Pivoting  15 

2.2.2  Jordan  Elimination  16 

2.3  Error  Analysis  21 

2.3.1  Basic  Operations  22 

2.3.2  Inner  Product  Computations  24 

2.3.3  Matrix  Operations  32 


CHAPTER  III  -  A  NORM  BOUND  ON  THE  ERROR  IN  COMPUTING 

A  REFLEXIVE  GENERALIZED  INVERSE 

3.1  Computation  of  a  Reflexive  Generalized 


Inverse  35 

3.2  An  Error  Analysis  40 

3.2.1  A  Bound  on  II F-,  II  42 

II  J_ll  OO 

3.2.2  A  Smaller  Bound  on  II  F^||  ^  52 

3.2.3  A  Bound  on  II  ^2 H  00  57 

3.2.4  A  Bound  on  II  -^3!!  oo  62 


3.2.5  A  Bound  for  the  Computed  Reflexive 
Generalized  Inverse 


65 


Page 


CHAPTER  IV  -  A  NORM  BOUND  ON  THE  ERROR  IN  COMPUTING 

THE  PSEUDOINVERSE 

4.1  Computation  of  the  Pseudoinverse  67 

4.2  An  Error  Analysis  70 

4.2.1  A  Bound  for  the  Computed  Hermite 

Normal  Form  70 

4.2.2  A  Bound  for  the  Product  P*AB*  77 

4.2.3  A  Bound  in  Computing  (P*AB*)”^  82 

4.2.4  A  Bound  for  the  Computed 

Pseudoinverse  87 

CHAPTER  V  -  NUMERICAL  RESULTS 

5.1  Test  Data  90 

5.1.1  Nonsingular  Matrices  90 

5.1.2  Rectangular  Matrices  with  Known 

Reflexive  Generalized  Inverses  94 

5.1.3  Rectangular  Matrices  with  Known 

Pseudoinverses  96 

5.2  Summary  of  the  Results  99 

5.2.1  Results  for  the  Reflexive 

Generalized  Inverse  100 

5.2.2  Results  for  the  Pseudoinverse  106 

CHAPTER  VI  -  CONCLUSIONS  AND  SUGGESTIONS  FOR 

FUTURE  RESEARCH 

6.1  Conclusions  112 

6.2  Suggestions  for  Future  Research  113 

BIBLIOGRAPHY  116 

APPENDIX  -  LISTINGS  OF  THE  FORTRAN  IV  SUBPROGRAMS  120 


LIST  OF  TABLES 


Page 


5.1 

Error  in 
Inverse 

the 

Reflexive  Generalized 

101 

5.2 

Error  in 

the 

Pseudoinverse 

107 

CHAPTER  I 


INTRODUCTION 


1 . 1  History 

The  concept  of  the  inverse  of  a  nonsingular  matrix  was 
first  generalized  to  include  all  matrices  by  Moore  [11]. 

The  generalized  inverse,  or  "general  reciprocal",  as  Moore 
called  it,  was  independently  rediscovered  by  Penrose  [13] 
in  1955.  Since  then,  many  persons  have  investigated  its 
properties.  Rado  [15]  showed  the  equivalence  of  the  defini 
tions  given  by  Moore  and  Penrose.  Other  definitions  have 
been  given  by  Rao  [16]  and  by  Goldman  and  Zelen  [4],  Some 
practical  methods  of  computation  for  the  different  types  of 
generalized  inverses  have  recently  been  suggested  by  Rao 
[16],  Pyle  [14],  Ben-Israel  and  Cohen  [1],  Graybill,  Meyer 
and  Painter  [5],  Rust,  Burrus  and  Schneeberger  [18],  and 
Willner  [23]. 

Following  Greville  [6]  and  Rohde  [17],  we  will  call 
the  generalized  inverse  as  defined  by  Moore  and  Penrose  a 
Pseudoinverse .  It  has  often  been  called  the  Moore-Penrose 
Generalized  Inverse  or  simply  the  Generalized  Inverse.  The 
Pseudoinverse  is  defined  by  the  following  theorem  due  to 
Penrose  [13].  We  denote  by  A*  the  conjugate  transpose  of 


a  matrix  A. 


rrari 
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Theorem  1.1  The  four  equations 


(1.1) 

AX A  =  A3 

(1.2) 

XAX  =  X3 

(1.3) 

(XA) *  =  XA 

and 

(1.4) 

(AX)  *  =  AX 

have  a  unique 

solution  for  any  matrix 

The  unique  solution  of  equations  (1.1),  (1.2),  (1.3), 
and  (1.4)  is  the  Pseudoinverse  of  A  and  is  written  as  A+. 

If  A  is  of  dimension  m-by-n,  then  A+  will  be  of  dimension 
n-by-m. 

Goldman  and  Zelen  [4]  have  proposed  the  following  general¬ 
ized  inverse. 


Definition  1 . 1  Let  A  be  any  real  matrix  of 


dimension  m 

-by-n.  Then  a  Weak 

Generalized  Inverse  of  A 

is  a  matrix 

A^n^  of  dimension 

n-by-m 

satisfying  equations 

(1.1)3  (1.2) 

and  ( 1 . 3) 3  i .  e  . 

(1.5) 

ii 

A, 

(1.6) 

A<n)  A  A(n) 

=  A(n) 

and 

(1.7) 

(A(n) A) '  = 

A(n)  A , 

3 


where  the  superscript  prime  in  (1,7)  denotes  matrix  trans¬ 
position. 

( n ) 

A  proof  of  the  existence  of  A  is  furnished  by 

Goldman  and  Zelen.  The  solution  of  (1.5),  (1.6)  and  (1.7), 

( n ) 

however,  is  not  unique.  The  notation  A  for  the  Weak 

Generalized  Inverse  is  due  to  Rohde  [17],  who  called  it 
the  Normalized  Generalized  Inverse. 

Rao  [ 1 6 ]  has  shown  the  existence  of  another  non-unique 
generalized  inverse,  for  which  Rohde  [17]  states  the  follow¬ 
ing  definition. 

Definition  1.2  A  solution  G  to  the  equations 

(1.8)  AG A  =  A 
and 

(1.9)  GAG  =  G 

is  called  a  Reflexive  Generalized  Inverse  of  a  matrix  A. 

Rao’s  proof  for  the  existence  of  a  Reflexive  Generalized 
Inverse  consists  of  a  constructive  procedure  which  may  be 
used  to  compute  a  Reflexive  Generalized  Inverse  (see  Chapter 
III)  o 

(  cr  ) 

A  Generalized  Inverse  A  &  of  an  arbitrary  matrix 
A  satisfying  only  equation  (1.1)  has  been  proposed  by  Rao 
[16].  Rohde  [17]  states  the  following  definition: 


X 
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(g) 


Definition  1.3  A  matrix  A 

(  n  ) 

A  =  A. 

to  be  a  matrix  such  that 


is  said  to  be  a 

(g) 


Generalized  Inverse  of  the  matrix  A  if  A  A 
Rao  originally  defined 
for  any  vector  y  for  which  Ax  =  y  is  consistent, 

(  cr  ) 

x  =  A  &  y  is  a  solution.  As  a  result  of  this  definition, 

(g) 


Rao  proves  that  A 
and  only  if  A  A^g^ 


is  a  Generalized  Inverse  of  A  if 

(g) 


A  =  A.  The  matrix  A 


is  also  not 


unique . 


1 . 2  Purpose  of  the  Study 

Wilkinson  [21]  has  stimulated  interest  in  the  study  of 
the  cumulative  effect  of  round-off  errors  when  performing 
computations  on  a  digital  computer.  An  important  criterion 
in  the  choice  of  a  suitable  method  for  the  computation  of  a 
particular  type  of  generalized  inverse  is  the  magnitude  of 
the  upper  bound  for  the  round-off  error  incurred  in  the 
computation.  In  this  paper  we  attempt  to  find  upper  bounds 
for  the  round-off  error  incurred  in  the  computation  of  the 
Pseudoinverse  and  a  Reflexive  Generalized  Inverse.  The 
computational  procedures  due  to,  Willner  [23]  and  Rao  [16], 
respectively,  are  examined. 


CHAPTER  II 


BASIC  THEORY 


2 . 1  Matrix  Algebra 

In  this  section  we  discuss  some  of  the  basic  properties 
of  matrices  which  are  used  throughout  this  paper.  Some 
particular  types  of  matrices,  which  appear  in  the  computation 
of  generalized  inverses,  are  defined.  The  concept  of  the 
norm  of  a  matrix  is  introduced,  and  various  types  of  matrix 
norms,  including  the  norm  of  a  rectangular  matrix,  are  also 
considered . 


2.1.1  Notation  and  Basic  Definitions  Upper-case  Latin 
letters  and  upper-case  Greek  letters  designate  matrices,  and 
all  lower-case  letters  designate  scalars.  We  say  that  a 
rectangular  matrix  A  =  (a„)  with  m  rows  and  n  columns 
is  of  dimension  m-by-n.  If  m  =  n,  the  matrix  A  is 
square  and  of  order  n.  The  zero  matrix  is  denoted  by  0 
and  the  identity  matrix  by  I.  The  complex  conjugate  of  a 
matrix  A  is  denoted  by  A*.  If  A  is  a  real  matrix.  A* 


is  written  as  A’ ,  the  transpose  of  A.  A 


-1 


designates 


the  inverse  of  a  nonsingular  matrix  A.  The  matrix  whose 
elements  are  the  absolute  values  of  the  corresponding  elements 
of  A  is  denoted  by  | A | ,  i.e.,  | A |  =  ( | a^ | ) . 


X  B  ’  B 
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For  a  square  matrix  A  of  order  n,  the  elements 
a^,  1=1, 2,..., n,  constitute  the  main  diagonal  of  A.  For 

a. ,  =  0  ,  i  j  , 
ij 

and  at  least  one  nonzero  element  of  the  main  diagonal,  A 
is  said  to  be  a  diagonal  matrix.  Moreover,  for 


•  -  o 

1  <  J 


5 


A  is  called  a  lower-triangular  matrix,  and  for 


1  >  j 


5 


an  upper-triangular  matrix. 

2.1.2  Permutation  Matrices  According  to  Wilkinson 
[22],  permutation  matrices  may  be  defined  as  follows: 

Definition  2.1  A  matrix  P  of  order 

————— —————— — — — — ^  1  J  0^  O  J  •  •  9  3  0^ 

1  2  n 

n  is  a  permutation  matrix  if 


a  . 


1 


and 


j  ^  ai  > 


• 
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where  (a  ,  ,  .  .  .  ,  )  is  some  permutation  of  ( 1 ,  2  ,  .  .  .  3n) . 

These  matrices  have  the  property  that  one  and  only  one 
element  of  every  row  and  column  is  nonzero  and  equal  to  1. 

The  permutation  matrices  used  in  this  paper  are  a 
subset  of  the  set  of  all  permutation  matrices  as  defined 
above.  They  are  obtained  by  interchanging  two  rows  of  the 
identity  matrix  of  order  n.  We  let  designate  the 

permutation  matrix  obtained  by  interchanging  rows  i  and 
J  of  I. 


Pre-multiplication  of  an  arbitrary  matrix  A  of 

dimension  m-by-n  by  a  permutation  matrix  R. .  has  the 

J 

effect  of  interchanging  rows  i  .  and  j  of  A.  Post¬ 
multiplication  of  an  arbitrary  matrix  A  of  dimension 

m-by-n  by  R. .  interchanges  columns  i  and  j  of  A. 

^  <1 

We  will  denote  a  permutation  matrix  used  for  post-multiplica¬ 


tion  by  C. .  rather  than  R. . . 
J  ij  ij 


2.1.3  Hermite  Normal  Form  Before  defining  the  Hermite 
normal  form  of  a  matrix,  we  will  define  what  is  meant  by  an 
elementary  row  operation. 


Definition  2.2  Three  elementary  row  operations 
may  be  defined  on  a  matrix  A  of  dimension  m-by-n: 

Type  I:  the  interchange  of  two  rows  of  A; 

Type  II:  the  addition  of  a  scalar  multiple  of  one  row 

of  A  to  another  row  of  A; 


- 
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Type  III:  the  multiplication  of  a  row  of  A  by  a 
nonzero  scalar. 

The  definition  of  the  Hermite  normal  form  of  a  matrix 
and  the  computational  procedure  are  due  to  Marcus  and  Mine 
[10]. 


Definition  2 . 3  Let  A  be  a  matrix  of  dimension 
m-by-n  with  rank  k.  Then  the  Hermite  normal  form  H  of 
A  is  a  matrix  of  dimension  m-by-n  such  that: 

(i)  the  first  k  rows  of  H  are  nonzero  and  the 
remaining  m  -  k  rows  are  zero ; 

(ii)  the  first  nonzero  element  in  the  i-th  row  of  H3 
i=l,23...3k3  is  1  and  occurs  in  column  n^} 
where  n  <  n ^  <  .  .  .  <  n and} 

(Hi)  the  only  nonzero  element  in  column  n ^  of  H  is 
the  1  in  the  i-th  row. 


Using  elementary  row  operations,  the  reduction  of  A 
to  H  can  be  performed  by  repeating  the  following  two 
steps  for  i=l,2,...,k. 


(i)  Let  n^  denote  the  column  number  of  the  i-th 
column  of  A  with  a  nonzero  entry  in  rows 
i,  i+l,...,m.  By  a  Type  I  row  operation  bring 
a  nonzero  entry  to  the  (i,n. )  position  from 
rows  i,i+l,...,m.  Then  make  this  element  equal 
to  1  by  a  Type  III  row  operation. 
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(ii)  With  Type  II  row  operations,  annihilate  all  of 
the  elements  of  column  n^  except  the  element 
(i ,n.) . 

The  computation  of  the  Hermite  normal  form  is  essential 
in  Winner’s  algorithm  for  the  computation  of  the  Pseudo¬ 
inverse.  We  have  programmed  the  above  procedure,  incorporating 
a  type  of  partial  pivoting  (see  Section  2.2). 

2.1.4  Matrix  Norms  A  norm  of  a  matrix  is  a  non¬ 
negative  scalar  which  is  a  measure  of  the  magnitude  of  the 
matrix.  A  norm  of  a  square  matrix  A,  denoted  by  ||  A  || ,  must 
satisfy  the  following  four  conditions: 


(i) 

II  A  II  > 

0  if  A 

/  0  and  ||  0 1|  =  0  ; 

(ii ) 

II  kA  ||  = 

1  k  |  II A II 

for  any  scalar  k  ; 

( iii ) 

II  A+B  || 

S  II  A  ||  + 

1 B  ||  ; 

(iv) 

II  ab  ||  < 

II  A||  ||  B  || 

• 

Pour  common  matrix  norms  defined  on  square  matrices  are: 


(2.1) 


max  l 

j  i 


max  l 

i  j 


j 


(2.2) 


5 
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(2.3)  ||a||2  =  (maximum  eigenvalue  of  A*A)^2  , 

and 

(2.M  |l  A||  =  (  l  l  laj-l2)172  . 

i  3  3 

Korganoff  and  Pavel-Parvu  [9]  show  that  these  definitions 
may  be  extended  to  Include  rectangular  matrices.  They  mention, 
however,  as  does  Householder  [8],  that  rectangular  matrices 

c 

can  be  normed  by  adjoining  zero  rows  or  columns  to  them  and 
then  applying  the  definitions  of  norms  of  square  matrices. 

This  is  how  we  shall  compute  the  norm  of  a  rectangular  matrix. 

2 . 2  The  Elimination  Methods  of  Gauss  and  Jordan 

There  are  various  well-known  elimination  procedures 
associated  with  the  name  of  Gauss  for  solving  a  system  of 
linear  equations  with  a  nonsingular  coefficient  matrix  and 
for  inverting  nonsingular  matrices.  In  this  section  we 
consider  two  of  these  elimination  methods,  Gaussian  elimina¬ 
tion  and  Jordan  (or  Gauss-Jordan )  elimination.  We  examine, 
in  particular,  the  matrix  equivalents  of  these  two  elimina¬ 
tion  methods  (see  Fox  [3]). 

2.2.1  Gaussian  Elimination  Let  A  be  a  matrix  of 
dimension  m-by-n,  where  m  <  n,  and  let  the  rank  of  A  be 
k,  k  <  m.  Then  A  can  be  reduced  to  a  matrix  of  the  follow¬ 
ing  form  (if  the  first  k  columns  of  A  are  linearly 
independent ) 


' 


■ 


11 


X  G 

' — *  r— I 

cd 


v — '  i — I 

cd 


cd 


X  cn 

v — "  i — I 

cd 


X  cm 

v — '  rH 

cd 


cd 


X  G 

CM 

cd 


w  CM 

cd 


cd 


J*  on 

^  CM 

cd 


X  CM 
^  CM 

cd 


X  G 
^  on 
cd 


^  on 
cd 


cd 


X  on 
^  on 
cd 


o 


^  G 
cd 


—  x 

cd 


X  * 

•  •  • 

o 

1 — 1 

^  CM 

^  on 

o 
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by  a  series  of  k  transformations. 

2. 2. 1.1  Partial  Pivoting;  By  partial  pivoting 
we  mean  the  following:  at  the  i-th  stage,  i=l,2,...,k. 


in  the  transformation  of  A  into  the  pivotal  element 


is  chosen  to  be  the  element  of  largest  magnitude  in  the  i-th 
column  and  in  rows  i,i+l,...,m.  If,  at  the  i-th  stage,  the 
pivotal  element  is  in  the  j-th  row,  then  rows  i  and  j 
are  interchanged.  This  interchange  can  be  accomplished  by 


R.  to  indicate  that  this  is  the  permutation  matrix  which 


i 


accomplishes  the  partial  pivoting  at  the  i-th  stage. 

In  terms  of  matrices,  the  partial  pivoting  at  the  first 
stage  may  be  represented  by  the  matrix  product  R^A.  Then 
the  first  reduced  matrix  will  be  J^R^A,  where  is  a 

lower-triangular  matrix  of  the  form 


1 


(2.6) 


J 


1 


0  1 


m 


ml 


0  0 


0  1 
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where 


(2.7) 


m 


il 


ll 

:11 


The  matrix  J^R  A  will  be  of  the  form 


(2.8) 


J1R1A  = 


a(1) 

a(1) 

a(1) 

a(D 

all 

12 

13 

*  •  CL 

In 

0 

a(1) 

22 

a(1) 

a(1) 

»  •  CL  /-\ 

2n 

0 

• 

a(1) 

a32 

• 

a(1) 

33 

• 

a(1) 

'  '  a3n 

• 

• 

0 

• 

a(1) 

am2 

• 

a(1) 

am3 

• 

a(1) 

mn 

(k) 

The  reduced  form  of  A,  A  }  can  be  obtained  from  the 
following  recursive  relations: 


(2.9) 


A(0)  =  A 


A(l)  =  J . R. A  ^ 1  1) 

l  l 


1=1,2 , . . . ,k, 


is  a  lower-triangular  matrix  of  the  form 


where  each 


in 


1 

0 

1 

0 

• 

0 

• 

• 

0 

• 

0 

.  .  .  1 

(2.10)  J.  = 

1 

0 

0 

0 

1 

0 

0 

.  .  .  0 

m.  ,  ,  . 

l  +  l ,  l 

1 

0 

• 

• 

0 

• 

• 

0 

• 

• 

i  +  2  ,i 

• 

• 

0 

• 

• 

• 

0 

• 

0 

• 

0 

• 

m  . 
mi 

e 

0  ...  1 

with 

(2.11) 

m .  . 
ji 

J1  t 

1=1,2 

j  =1+1 

5  .  .  .  j  k  j 

,i+2  j . . . ,m, 

a(i-D  ’ 
11 

where  af ?  ^  = 

air 

From  (2.9), 

the  final 

reduced  matrix 

(2.12) 


J,  R,  J,  n  R, 
k  k  k-1  k-1 


*^1^1^  * 


(k ) 

This  process  of  reducing  A  to  A  is  called  Gaussian 

forward  elimination  with  partial  pivoting.  As  a  result  of 


.  . 


15 


the  partial  pivoting,  we  note  that 

1=1*2, . . . ,k, 

J =i+l ,i+2 , . . . ,m. 


(2.13) 


m 


ji 


<  1 


2. 2.1. 2  Complete  Pivoting  An  alternative  to 
partial  pivoting  is  complete  pivoting,  which  may  be  described 
as  follows:  at  the  i-th  stage,  i=l,2,...,k,  in  the  trans¬ 
formation  of  A  into  A^\  the  pivotal  element  is  chosen 
to  be  the  element  of  largest  magnitude  in  columns  i,i+l,...,n 
and  in  rows  i,i+l,...,m.  In  order  to  transfer  the  pivotal 
element  to  position  (i,i)  of  the  matrix,  one  row  interchange 
and  one  column  interchange  are  necessary.  If  the  pivotal 
element  at  the  i-th  stage  is  a^1-  ,  then  rows  i  and  r 
and  columns  i  and  s  must  be  interchanged.  This  can  be 
accomplished  by  pre-multiplication  by  R^  and  post-multi- 

plication  by  C.  .  We  shall  omit  the  double  subscript  and 

1  s 

simply  write  and  C  . 

Thus  the  complete  pivoting  at  the  first  stage  can  be 
represented  by  the  matrix  product  R^AC^.  Then  the  first 
reduced  matrix  will  be  J  R^AC^,  where  is  defined  by 

(2.6),  and  will  be  of  the  same  form  as  J^R-^A.  The  final 
reduced  form  of  A  can  be  obtained  from  the  recursive 


relations 


. 
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(2.14) 


A(0)  =  A 


A(l)  =  J.R.A(l  1)C.J  i=l , 2 , . . . ,k, 

11  l 5  9  9  99 


where  the  matrices  J  are  defined  by  (2.10)  and 
a ^  .  The  final  reduced  matrix  is 


(2.15) 


A 


(k)  =  Wk-lRk-r--JlRlAClC2---Ck  • 


(k) 

This  process  of  reducing  A  to  A  is  called  Gaussian 

forward  elimination  with  complete  pivoting. 

As  was  the  case  for  partial  pivoting,  it  may  be  noted 

that 


(2.16)  Irru^l  <  1  ,  i=l,2,...,k, 

j =i+l,i+2 , . . . ,m. 

With  complete  pivoting,  the  following  inequality  also  holds 
on  the  elements  of  : 


(2.17) 


l  1 , 2  ,  .  .  .  ,  , 

j  =i+l ,i+2 , . . . ,n. 


2.2.2  Jordan  Elimination  We  again  let  A  be  a  matrix 
of  dimension  m-by-n  with  rank  k,  where  m  <  n  and  k  <  m. 
It  is  well-known  that  Jordan  elimination  reduces  a  nonsingular 


17 


matrix  of  order  n  to  diagonal  form  in  n  transformations 
(see  Fox  [3]).  When  applied  to  a  rectangular  matrix  A,  k 
Jordan  transformations  reduce  A  to  A^k\  where 


(2.18)  A 


(k) 


(k) 

'll 


0 


0 


0 


(k) 

'22 


0 


0 


0 


kk 


(k) 

*l,k+l 


(k) 

'2  ,k+l 


(k)  (k) 


k  ,k+l 


(k) 

'In 


(k) 

'2n 


‘8° 

kn 


0 


if  the  first  k  columns  of  A  are  linearly  independent. 

The  first  Jordan  transformation  is  identical  to  the 
first  transformation  of  Gaussian  elimination.  In  terms  of 
matrices,  it  can  be  represented  as  the  matrix  product  J^A, 
where  J^  is  defined  by  (2.6).  All  succeeding  Jordan 
transformations  differ  from  the  corresponding  Gaussian 
transformations.  The  second  Jordan  transformation  can  be 
represented  by  the  matrix  product  J^J-^A,  where 
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(2.19) 

with 

(2.20) 

and  A 

(2.21) 


J2  = 


m 


12 


0 


0  m 


0  m 


32 


42 


0  m 


m2 


0  0 


0  0 


0  0 


0 


0 


1  0  .  .  .  0 


0  1  .  .  .  0 


m12 


a(1) 

12 

Tnr  * 

a22 


1  =  1 , 3 , 4 , . . . ,m. 


=  J^A.  Then  will  be  of  the  form 


J2J1A  = 


a 


(2) 

11 


0 


0 


0 


0 


'22 


0 


0 


(2) 

'13 


(2)  (2) 


■23 


(2) 

'33 


l(2) 

m3 


(2) 

'In 


(2) 

'2n 


(2) 

'3n 


t(2) 

mn 


(k) 


can  be  obtained  from  the  following  recursive 


relations : 


' 
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(2.22) 


A<°>  -  A 


A 


(l)  =  JiA(l“1)  , 


1=1 ,2 , . . .  ,k 


where  each  matrix  J  has  the  form 


(2.23)  J .  = 


l 


1  0  ...  0  m 


0  0 


0  0 


11 


0  1  ...  0  m 
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0  ...  1  m.  . 

i-l5i 


0  0  ...  0  1 


0  m 


1+1,1 


0  m 


ml 


0  ...  0 


0  ...  0 


0  ...  0 


0  ...  0 


0 


0 


with 


(2.24) 


m .  .  = 
Ji 


- ( i-l) 

cL  «  « 

■ii 

ital  > 

CL  •  • 

11 


1=1,2, . . . ,k, 

j  =  1 , 2 ,  .  .  .  ,  m , 


and 


The  final  reduced  matrix  is  given  by 
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(2.25) 


A 


(k) 


J  T  t  t  c  •  e  ti  . 

k  k-1  1 


In  practice,  some  form  of  pivoting  may  be  necessary 

to  insure  that  none  of  the  pivotal  elements  af^-^  become 

n 

zero.  However,  to  achieve  the  form  of  A in  (2.18),  the 
pivot  at  the  i-th  stage  must  be  chosen  from  rows  i , i+1 , . . . ,m . 
Thus  it  follows  that 

(2.26)  |m  |  <  1  ,  i=l , 2 , . . . ,k , 

J 

j =i+l ,i+2 , . . . ,m, 

but  the  magnitudes  of  the  multipliers  m.„,  j =1 , 2 , . . . , i-1 , 

J  ^ 

are  not  bounded. 

By  analogy  with  the  method  of  Gaussian  elimination,  the 
matrix  equivalent  of  Jordan  elimination  with  partial  pivoting 
is 


(2.27) 


J,  R,  J,  R,  -)...J-IR-IA  , 
k  k  k-1  k-1  1  1  5 


and  the  matrix  equivalent  of  Jordan  elimination  with  complete 
pivoting  is 


J,  R,  J,  R, 

k  k  k-1  k-1 


. . j1r1ac1c2 


C 


k  ’ 


(2.28) 


•  •  • 
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In  both  (2.27)  and  (2.28),  however, 
on  the  magnitudes  of  the  multipliers 
the  main  diagonal. 


we  cannot  place  a  bound 

,  m . . ,  of  J .  above 
5  0 1 5  i 


2 . 3  Error  Analysis 

In  this  section  we  state  some  of  the  fundamental  floating¬ 
point  bounds  on  the  round-off  error  incurred  in  performing 
certain  mathematical  computations  on  a  digital  computer.  Only 
floating-point  computation  is  considered,  as  fixed-point  compu¬ 
tation  is  not  suitable  for  storing  numbers  of  a  sufficiently 
large  range  of  magnitudes  or  for  performing  lengthy  calcula¬ 
tions  with  numbers  of  varying  magnitudes.  All  of  the  results 
presented  here  are  due  to  Wilkinson  [21]. 

The  number  of  binary  digits  in  a,  the  mantissa  of  a 
floating-point  number  x  =  2  (a),  is  denoted  by  t.  Through¬ 
out  this  paper  we  assume  that  the  number  e  is  the  unit 
round-off  error,  i.e. 

(2.29)  | e |  <  2_t  . 

We  will  be  concerned  only  with  the  so-called  "backward 
error  analysis".  In  performing  a  backward  error  analysis  on 
a  certain  computation,  the  actual  difference  between  a  computed 
value  and  a  true  value  is  not  sought.  Instead,  where  our 
mathematical  computation  takes  the  form  x  =  f (a^ ,a2 , •  •  • >an) , 
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we  attempt  to  show  that  the  computed  value  obtained  for  x 
is  exactly  equal  to  f  (a^+e-^  }a^+e^  , .  .  .  >an+en)  for  some 
values  of  the  e^.  The  problem,  then,  is  to  place  bounds 
on  the  magnitudes  of  the  e^.  The  mathematical  equation 
which  defines  the  computed  value  of  x  is  called  a  compu¬ 
tational  equation,  and  we  write 

(2.30)  x  =  f (a,+e, ,a«+e0, . . . ,a  +e  )  , 

1  1  ’  2  2 5  5  n  n  5 

followed  by  inequalities  which  limit  the  magnitudes  of  the 

e .  . 

1 

bi 

2.3.1  Basic  Operations  Let  x  =  2  (a.)  and 

b2 

x^  =  2  (ap  be  two  floating-point  numbers,  and  let  fl 
denote  a  floating-point  arithmetic  computation.  Then  the 
following  theorem  gives  a  bound  to  the  round-off  error 
incurred  in  performing  a  floating-point  addition,  subtraction, 
multiplication  or  division  with  a  double-precision  accumulator. 

Theorem  2.1 

fl(x^Xx^)  =  (x^Xx^)(l+z)  , 

where  X  is  any  one  of  the  four  arithmetic  operators  +3 
x,  or  v. 
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Proof :  Let  the  exact  normalized  value  of  x-^Xx^ 

b3 

be  2  3 (a^ ) .  Then 

flCx^Xx^)  =  (x^Xx2)  +  e1  , 

where 

|e1|  <  |  2_t(2b3)  . 

Since  a  relative-error  bound  of  the  form 


fl(x1Xx2)  =  (x1Xx2) (l+e2)  , 


is  required,  we  can  write 


(x1Xx2)  e2  =  e1  . 


Therefore 


1 


x^Xx2 


\  2-t ( 2  3) 

|  (2b3) 


as 


due  to  normalization  of  floating-point  numbers. 
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Thus 


Q .  E  .  D . 


2.3.2  Inner  Product  Computations  Whenever  it  is 

necessary  to  compute  an  inner  product,  we  will  assume  that 

it  is  done  by  using  floating-point  accumulation,  which  may 

be  defined  as  follows.  If  the  inner  product  to  be  calculated 
n 

is  £  a.b.,  then  each  product  a.b.  is  computed  by  using 

i  =  l 

single-precision  floating-point  arithmetic  and  storing  the 
2t-digit  product  without  rounding.  Then  the  n  2t-digit 
products  are  added  using  double-precision  arithmetic,  and 
the  final  result  is  rounded  to  t  digits .  We  shall 
indicate  the  computation  of  an  inner  product  using  floating¬ 
point  accumulation  by  the  symbol  fl^.  The  symbol  fl^ 
denotes  a  floating-point  arithmetic  computation  using  a 
single-precision  accumulator. 

The  following  lemmas  are  required  for  the  proof  of 
the  next  theorem  (Theorem  2.2). 

Lemma  2 . 1  Using  a  single-precision  (t-digit) 
accumulator , 


=  (x i+x  g) ( l+£  j)  3 


f  l  ^  ^  ^  2~^x  2  ^ 


- 
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where 


e  I  <  —  2~^ 
1  1  ~  2 


Proof :  When  using  a  single-precision  accumulator, 

the  maximum  amount  of  round-off  error  will  be  incurred  when 
two  separate  rounding  errors  are  made  in  performing  the 
addition.  One  rounding  can  occur  before  the  addition  is 
carried  out,  in  shifting  the  smaller  of  the  two  floating¬ 
point  numbers  to  the  right  to  align  the  binary  points.  We 
will  denote  this  rounding  error  by  n^.  A  second  rounding 
error  may  be  incurred  after  the  addition  takes  place,  if  the 
sum  of  the  mantissas  is  greater  than  one.  In  this  case,  the 
mantissa  of  the  sum  must  be  normalized;  we  will  denote  this 
rounding  error  by  rip. 

bl  b2 

Let  the  exact  sum  of  x^  =  2  (a^)  and  x^  =  2  (a^) 

b  t 

be  2  (a^)j  and  assume  that  x^  <  Xp.  Then 


1  -t 

n-J  *  j  2  b ( 2  b  ) 


and 


n2l  s  |  2  t(2  3) 


Therefore , 
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n1+n2l  ^  |n1l  +  |n2 


2  -t  t’o-! 

<  |  2  t(2  J  ) 


S  I  2_t  ! xl+x2 I  5 


as  the  round-off  error  n2  will  not  arise  unless 

i  i  1  ^3 

|x^+x2|  >  ^  2  .  Since 


fi^(x1+x2)  =  (x1+x2)  +  n1+n2  j 


if  we  write 


fl^(xl  +  x2)  =  ( X 1  +  X  2 )  ( 1  +  £ 1  )  , 


then 


Q .  E .  D . 


Lemma  2 „ 2  Using  a  double-precision  accumulator 
but  without  floating-point  accumulation 3 

n  n 

fl(  \  a.b.)  =  £  a.b.(l+E.)  > 

7s  7s  ^  7s  7s  7s 
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where 


(l-2't)n  <  l+z1  <  (1+2  t)n 


and 


(1_2-t)n-r+2  s  1+£ 


<  (i+2-t)n-r+Z 


r=2 3 3  j .  .  . j  n . 


Proof:  We  define  quantities  s^  and 


recursively 


by  the  relations 


p  =  fl(a  b  ) 


and 


sn  =  p,  .  s  =  f  1  ( s  n  +p  )  . 
1  ^l5  r  r-1  ^r 


Then  it  follows  from  Theorem  2.1  that 


a  b  (l  +  £  )  , 
r  r  r  5 


where 


< 


2 


-t 


5 
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and 


sr  =  (sr-l+pr)(1+r,r)  ’ 


where 


Hence 


where 


n 


T  a. b . ( 1+e .  ) 
s*  11  1 


3 


1+e-L  =  (l+C1)(l+n2)  •••  (l+nn) 


and 


l+er  =  (l+^r ) (l+nr )  ...  (l+nn),  r=2 , 3 


3 


n. 


Thus 


(l-2"t)n  <  l+e1  <  (l+2“t)n 


and 


(l_2-t)n-r+2  s  1+£^  s  (1+2-t}n-r+2 


r  2,3,. . . ,n . 


Q .  E .  D  . 
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Theorem  2„2 


(2.35) 


n 


fl2(  I  a^.) 


n 

y  a  .b  .)  ( 1+z) 


n 

I 


a  .b  .e  . 

t  t  ^ 


3 


where 


(2.36) 


r 


<  |e  I  <  4  (n-r+2)  2  ^^2 
'  r  1  2 


Proof :  It  follows  Immediately  from  Lemmas  2.1  and 


20 2  that 


( 2  „  37  ) 


n 


f 10 (  7  a. b .  )  = 

2  ^  li 


n 

I 


aibi (l+ei ) 


( 1+e )  5 


where 


(2.38) 


(1-|  2’2t)n  s  1+e  <  (1+|  2"2t)n  . 


and 


(2.39) 


(1-|  2-2t}n-r+2  £  1+£^  s  (1+|  2-2t}n-r+2 


r=2,3. 


Jn 


O 


•  •  • 


t 
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In  the  result  of  Lemma  2.2,  the  term  2-t  is  replaced 
_  2t 

by  2  as  the  products  a^b^  are  stored  as  double-precision 

numbers.  The  factor  — ■  is  introduced  since  n  double¬ 
precision  numbers  a^b^  are  added  using  double-precision 
arithmetic  (replace  "single-precision”  by  "double-precision" 
in  Lemma  2.1).  The  factor  (1+e)  accounts  for  the  rounding 
of  the  double-precision  sum  to  single-precision. 

If  we  restrict  a  number  p  such  that 

|  p  2-2t  <  0.1  , 

then  expansion  of  (1  ±  2”  ;p  according  to  the  binomial 

theorem  leads  to  the  results 

(1+|  2-2t)P  <  1+|  p ( 1 . 06  )  2-2fc 

and 

(1-|  2~2t)p  >  1-|  p ( 1 . 06 )  2_2t  . 

If  we  define  t^  by 


(1.06)  2“2t  =  2"2t2  , 


then  (2.38)  and  (2.39)  may  be  written  as 
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(2,40) 


e  I  <  1  n2“2t2 

l1  ’  2 


and 


(2.41) 


e  |  <  |  (n-r+2)  2~2t2 
r 1  2 


respectively,  assuming  that 


|  n2  2t  <  0.1 


(2.37)  may  be  written 


n 


n 


n 


(2.42)  flp(  l  a.b. )  -  (  l  a.b.)(l+e)  =  (  l  a. b . e „ ) (1+e ) . 

^  1  1  T  .  11  T  111 


The  factor  1.06  introduced  in  the  definition  of  t^  is 
sufficiently  generous  that  the  factor  (1+e)  on  the  right- 
hand  side  of  (2.42)  may  be  dropped.  Thus 


n 


n 


(2.43) 


fl2(  I  a.b.) 


(  l  a.b. )(l+e) 
1 


n 

y  a. b . e  .  . 

“  ill 


Q .  E .  D . 


Since  the  in  Theorem  2.2  are  of  the  order  e  ,  they 

may  be  neglected,  for  a  sufficiently  large  value  of 
(2.43)  becomes 


t.  Then 
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n  n 

(2.44)  fl2(  l  a.b.)  e  (  l  a.b1)(l+e)  . 

_  2 1 

The  neglecting  of  error  terms  of  the  order  2  does  not 

usually  affect  the  over-all  result  significantly,  since  t 
is  greater  than  20  for  most  large  computers  today. 

2.3.3  Matrix  Operations  The  final  bound  we  wish  to 
mention  is  a  bound  on  the  round-off  error  incurred  in  a 
matrix  multiplication.  As  each  element  of  a  matrix  product 
is  an  inner  product,  equation  (2.44)  may  be  used  to  prove 
the  following  theorem. 

Theorem  2 . 3  If  A  and  B  are  matrices  conform¬ 
able  for  multiplication ,  then 

(2.45)  fl2(AB)  E  AB  +  E  , 

where 

(2.46)  |ff|  <  e  \A\  | B | 

Proof :  Let  C  =  AB.  Then 

n 

cij  =  k|1aikbkj  ’ 

where  n  is  the  common  dimension  of  A  and  B,  and  by  (2.44), 
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fl0(c. . 
2  1 3 


)  e 


n 

(  I 

k=l 


aikbkj  + 


Hence 

fl2(C)  =  fl2(AB) 

=  (1+e)  AB 
=  AB  +  eAB  . 

Therefore,  a  comparison  of  the  above  equation  with  (2.45) 
gives 

|  E  |  <e  |  A  |  |  B  | 

Q .  E .  D  . 


As  a  norm  bound,  (2.46)  can  be  written 


(2.47)  II E  II  s  £  II  A  III  B  II  . 

Whenever  a  multiplication  or  division  is  performed  using 

powers  of  2  (in  a  binary  machine),  no  round-off  error  is 

incurred.  This  is  true  as  only  the  exponent  need  be  changed, 

and  not  the  mantissa.  For  example,  in  the  computation 
a 

y  =  x  X  2  ,  where  X  represents  either  of  the  floating-point 


operations  of  multiplication  or  division,  the  mantissa  of  y 
equals  the  mantissa  of  x,  and  the  exponent  of  y  is  computed 
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from  a  and  the  exponent  of  x.  Consequently,  a  matrix 
multiplication  involving  any  permutation  matrix  does  not 
produce  any  additional  round-off  error. 


CHAPTER  III 


A  NORM  BOUND  ON  THE  ERROR  IN  COMPUTING  A 
REFLEXIVE  GENERALIZED  INVERSE 


3. 1  Computation  of  a  Reflexive  Generalized  Inverse 
An  elimination  method  for  computing  a  Reflexive 
Generalized  Inverse  G  of  an  arbitrary  m-by-n  matrix 
is  given  by  Rao  [16].  We  will  assume  that 

(3.1)  m  <  n  . 

If  m  >  n,  we  will  consider  the  transpose  of  A,  since  it 
can  easily  be  shown  that  a  generalized  inverse  of  A’  is  G’. 
A  procedure  for  computing  G  is  outlined  below. 

Partition  A  such  that 

(3.2)  A  =  [L | R]  , 


where  L  is  a  square  matrix  of  order  m  and  R  is  of 
dimension  m-by-(n-m).  Then  reduce  L  to  an  upper-triangular 
matrix  by  applying  Gaussian  forward  elimination  with  complete 
pivoting  to  the  rows  of  A,  obtaining 
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(3.3)  A 


(1) 


where  the  elements  x  represent  components  which  are  not 
necessarily  zero  and  where  r,  r  <  m,  is  the  rank  of  A. 
Next  apply  Gaussian  forward  elimination  to  the  columns  of 
A^1^  to  obtain  a  matrix  of  the  form 


(3.4) 


dl 

0 

0 

o 

• 

o 

© 

0 

0 

d2 

0 

0 

0 

0 

0 

on 

rd 

0  •  ©  o 

0 

0 

• 

• 

0 

• 

c 

• 

0 

• 

• 

• 

0 

• 

0 

o 

• 

0 

0 

0 

0 

d 

r 

0 


In  practice,  the  forward  elimination  process  which  transforms 
A^1^  into  A  need  not  be  carried  out,  since  we  can  simply 


write  (3.4)  from  (3.3). 
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Let  us  define  A  ,  a  matrix  of  dimension  m-by-na  as 


Let  us  also  define  square  matrices  ,  i=l,2,...,r,  of 
order  m  as  follows:  each  matrix  P^  is  identical  to  the 
unit  matrix  except  for  the  elements  in  its  i-th  column  below 
the  main  diagonal.  These  elements  are  the  multipliers  from 
the  Gaussian  forward  elimination  which  transforms  A  into 
For  example, 


(3.6) 


1  0  0  0  .  .  .  0 

0  1  0  0  .  .  .  0 


0  m^2  1  0 
0  m^2  0  1 


0 

0 


0 


0 


0 


1 


For  r  =  m,  P  =  I. 


38 


Similarly,  let  us  define  square  matrices  ,  i=l,2,...,r, 
of  order  n,  such  that  each  matrix  is  identical  to 

the  unit  matrix  except  for  the  elements  in  the  i-th  row  to 
the  right  of  the  main  diagonal.  These  elements  are  the 
multipliers  required  for  the  Gaussian  forward  elimination 
in  transforming  A ^ ^  into  A.  For  example, 

1  0  0  0  ...  0 

0  1  m  ^  ^  ^  ||  •  •  •  m  ^  ^ 

0  0  1  0  .  •  .  0 


(3.7)  Q2  = 


0  0  0  0  ...  1 

For  m=n=r,  Q  =  I. 

Let  R  and  CL,  i=l,2,...,r,  be  permutation  matrices 

of  orders  m  and  n,  respectively,  which  accomplish  the 
complete  pivoting  at  the  i-th  stage  of  the  forward  elimina¬ 
tion  in  reducing  A  to  A^1^.  Let  us  define 


(3.8) 


P  =  P  R  P  n  R  .  .P-.R-, 
r  r  r-1  r-1  1  1 


and 


r 


(3.9) 


Q  —  C C 2  •  •  *  ^pQ]_Q2 


.  .  .  Q 
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Then  P  and  Q  are  non-singular  matrices  such  that 


(3.10) 


PAQ  =  A. 


Theorem  3 . 1  A  Reflexive  Generalized  Inverse  of 
A  is  given  by 

(3.  ID  G  =  Q( A“;  'P  . 


Proof:  By  (3.10), 


(3.12)  A  =  P"1AQ”1  . 

Therefore , 

AGA  =  P"1AQ_1Q(A'")  ,PP"1AQ"1 
=  P~1A(A")’AQ"1 
=  P“1AQ"1 
=  A. 

Also , 
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GAG  =  Q(A“) ,PP"1AQ”1Q(A") ’P 
=  Q(A“) ’ A(A“) ’P 
=  Q(A“)’P 

G*  Q.E.D. 


3 . 2  An  Error  Analysis 

Let  F^,  Fp  and  F~,  be  the  matrices  of  round-off 
errors  incurred  in  the  calculation  of  Q,  (A-)’  and  P, 
respectively.  Then  F^  and  F^  are  square  and  of  orders 
n  and  m,  respectively,  and  F^  is  of  dimension  n-by-m. 

The  computational  equation  for  the  calculation  of  G  becomes 


(3.13) 


G 


(Q+F1)((A”),+F2)+F 


(P+F3)+F5 


or 


(3.14)  G  e  Q(A")*P+ 


F1  ( (A“)  ’+F2)+QF2 


P+(Q+F1)((A’)i+F2)F3 


+F4(P+F3)+F5  , 


where  F and  Fc  are  the  matrices  of  round-off  errors  due 
to  the  multiplications  of  Q  by  (A-)’  and  Q( A  )'  by  P, 
respectively.  We  will  denote  the  exact  value  of  G  by  G£ 


' 
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and  the  computed  value  of  G  by  Gc.  Then  it  follows  from 
(3.14)  that 


(3.15)  |  G^— G^, 


FX((A  )»+F2)+QF 


P+(Q+F1) ( (A”) ’+F2)F3 


+f4(p+f3)+f5 


Thus  we  obtain  the  following  norm  bound  on  the  error  incurred 
in  the  calculation  of  G: 


(3.16) 


We  now  proceed  to  place  bounds  on  ||  F  ||  ,  1  =  1,253,4,5. 
The  following  norm  will  be  used,  since  it  can  easily  be 
computed : 

(3.17)  II  A  Hcc.  =  max  l  I  ai  1  I  » 

i  j  J 

where  A  =  (a^ )  is  an  arbitrary  square  matrix.  Equation 
(3.17)  can  be  extended  to  include  rectangular  matrices,  since 
these  can  be  made  square  by  bordering  them  with  zeros  (see 
Section  2.1). 


. 
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3.2.1  A  Bound  on  ||  P  ||  ^ 
computed  values  of  Q  by  Q„ 
Therefore 

(3.18)  I qe_QC 1 


Let  us  denote  the  exact  and 
and  ,  respectively. 


Let 

(3.19)  Q(1)  =  Q.,Q2.  .  .Qr  . 

Then  the  calculation  of  Q„  from  qU)  Is  accomplished  by 
appropriate  row  interchanges,  and  involves  no  arithmetic 
computations.  Thus  if  we  denote  the  exact  and  computed 
values  of  by  and  Q^,^,  respectively,  it 

follows  that 

(3.20)  IIQe1)-QC1)1o°  =  1!  FlH  oo  • 


In  order  to  place  an  upper  bound  on  the  round-off  error 
incurred  in  the  calculation  of  let  us  define  recur¬ 

sively  matrices  T^  as  follows: 


(3.21) 


T1  =  Q 

T^  =  j_  +  ]_^  ]_j  i  —  2,3j...ji1. 

i 


■ 
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Then 


(3.22) 


First  an  upper  bound  on  the  error  incurred  in  the 
matrix  multiplication  Q  .  ,,T.  ,  must  be  obtained.  All 
of  the  elements  of  the  matrices  ,  i=l,2,...,r,  are 

less  than  or  equal  to  one  in  magnitude  since  complete 
pivoting  is  used  to  transform  A  into  A*'1^.  Hence 


(3.23)  | Q± |  * 


1  0...0  0  0  0...0 


0 


0  0 


0  0 


0  0 


0  0  0  0  .  .  .  0 


0 


0  0  0 


0  0  ...  0  0 


1 


0 


0  0  0  0 


0 


1 


0 


where  only  the  i-th  row  of  the  matrix  in  (3.23)  contains 
more  than  one  nonzero  element.  Using  (3.23),  it  can  easily 

be  shown  that  the  product  I ^r-i  +  1 ^ ^r-i  +  2 I ' *  *  I I  ls  bounded 
by  the  matrix  given  in  (3.24),  and  hence,  | T± |  is  also  so 


bounded . 
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In  the  matrix  (3.24),  only  rows  r-i+1  to  r,  inclusively, 
contain  more  than  one  nonzero  element. 

Let 


0 


•r-i+1 


=  (qij) 


=  <V 


) 


and 


Qr-i+lTi-l  =  Ti  “  > 


From  the  structures  of  these  matrices,  it  can  be  seen  that 


(3.25) 
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Thus,  in  the  multiplication  of  1  by  Qr_i+i,  only  the 

elements  in  row  r-i+1  of  to  the  right  of  the  main 

diagonal  need  be  computed,  the  remainder  being  identical 
to  the  corresponding  elements  of  Therefore,  the 

matrix  ^  of  round-off  errors  for  the  multiplication 

of  T^  i  by  is  null  except  for  row  r-i+1.  Let 


0 


0 


0 


0 


0 


(3.26)  E._1 


er-i+l,r-i+2  er-i+l,r-i+3  *'*  er-i+l,n 


0 


0 


0 


0  .  .  .  0 


, .  .  .  ,  n 


i  J 


is  computed  as  an  inner  product  of  a  row  vector  of  Q  .+1 
and  a  column  vector  of  T^  We  now  use  the  result  of  (2.44) 

to  place  a  bound  on  the  nonzero  elements  of  (3.26).  From  (3.24) 
it  can  be  seen  that 


(3.27) 


<  2 


i-1 


j  =r-i  +  2 ,r-i  +  3 


r-i+1 J 


Therefore,  we  have 


<  2 


i-1 


n . 


e 


r-i+1, j 


(3.28) 


j  =r-i  +  2 ,r-i+3 , 


•  •  • 


- 
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The  following  theorem  is  a  result  of  equation  (3.28). 

Theorem  3*2  Let  Q  .  3  i=l32t...3r3  be  square 

'  'Is 

(2  ) 

matrices  of  order  n  which  transform  A  into  A.  Define 

matrices  T ^3  i=l323...3r3  by  (3.21).  Then  for 


(3.29) 


T  . 


Q  ,  •  ,  7  T  .  7  +  E  .  n  3 
r+v+1  i-l  i-l 


we  can  bound  the  matrix  E.  of  round-off  errors  as  follows 

Is  —  1 


(3.30)  ||  S'  •  _  7 II  oo  *  (n-r+i-1) 2Z~1 z3  i=23  33  ...  3  r 


i-l 


is  calculated  recursively  as  follows 


(3.3D  1  . 


T  =  Q 
1  r 


T  n  =  Q  Q  +  E , 
2  r- 1  r  1 


Q  oQ  i Q  +  Q  0E,  +  E0 
r-2  r-1  r  r-2  1  2 


T^  -  Q-^Qp  •  •  •  t  Q-|  Qp  •  •  •  pE-|  +  Q -|  Q p  »  *  *Q-|o_pE 


+  ... 


1^2  *  *  *  'r-2  1  12  *  •  *  r-3  2 


+  QiQ2Er-3  +  QlEr-2  +  Er-1 


and 


. 


' 


. 
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(3.32) 


Therefore 


(3.33)  IQe1)"QC1)I  =  lQiQ2‘  ’  *Qr>-2Ei  +  QnQ?- •  .Qr_^EP+.  .  . 


r-2  1 


12 


‘r-3  2 


+  Q, Q0E  0  +  Q  E  0  +  E  , 
1  2  r-3  1  r-2  r-1 


Using  (3.23),  it  may  be  shown  that 


< 

. 
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Hence , 


(3.35) 


llQ1IUlQpL---«QiIL  s  n2 


i-1 


i  =  l  ,2 , . . . ,r 


Using  (3.33), 

(3.36)  IIQ^-Q^IL  *  iQjjQpIL-- JWIJEJL 

+  II  Qj  Jq2IL-  •  4<ar_3ll  Je2IL+-  •  --HIqJ  jQ2IIJEr_3 
+  ll«iUV2l-+lEr-ll.  • 

Using  (3.30)  and  (3.35)  in  (3.36), 

(3.37)  llQE1)_QC1)Noo  -  2r~3n(n-r+l)2e 

+  2r“^n (n-r+2 ) 4e+ . . . +2n (n-3 ) 2r  3e 
+  n(n-2)2r-2e  +  (n-l)2r_1e 

=  n2r_2e  £(n-r+l)  +  (n-r+2) 

+  ...+  (n-3)  +  (n-2)  +  2(^I 

<  2r"’3n(r-l)  (  2n-r )  e  , 


■ 
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where  the  last  Inequality  holds  only  if  n  >  2.  The  following 
theorem  is  a  result  of  equations  (3.20)  and  (3.37). 

Theorem  3 . 3  Let  F ^  be  the  matrix  of  round-off 
errors  in  the  calculation  of  Q  according  to  equation  (3.9). 
Then 


(3.38) 


||  Fjoo  <  n(r-l  )  (  2n-r)  z 


for  n  >  2. 


3.2.2  A  Smaller  Bound  on 


Ft  II 


1"  00 


We  now  attempt  to 


place  a  smaller  norm  bound  on  by  bounding  the  right-hand 
side  of  equation  (3.33),  using  (3.26)  and  (3.3*0.  From  (2.44) 
and  (3.24),  it  follows  that  (3.26)  may  be  bounded  by 


(3.39) 


l-l 


<  £ 


0 


0 


0  0  0  0 


0  ...  0 


0 


0 


0 


0 


0 


,i-2 


0 


0 


,1-1 


0 


0 


,i-l 


0 


0 


,i- 


0 


— . . — 
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where  only  the  elements  In  the  (r-i+l)-th  row  and  columns 
r-i+2,  r-i+3s...,n  are  nonzero.  Direct  multiplication  of 
(3.3*0  and  (3.39)  produces  the  following  result: 


(3.40) 


Q1 I  I Q2 


Qi I ' Er-i- 1 


<  e 


0  ...  0  21"1  21 


0 


0 


0  .  .  .  0 


0  .  .  .  0 


0  .  .  .  0 


0  .  .  .  0 


2i+1  ...  2r~2 


0  .  .  .  0  21  2  21”1  21  ...  2r“3 


1 


1 


0 


0 


,r-i 


02  9r-i-l 

C—  •  o  •  C— 


22  ...  2r-1"1 


0  0 


0  0 


0 


0 


2r“2  ...  2r"2 


,r-3 


.  .  .  2 


r-3 


2r-i  ...  2r~1 


2r-1'1  ...  2r"i"1 


r— i— 1  r— i— 1 

2r  ^  1  ...  2r  1  1 


0 


0 


0 


0 


where  rows  i+2 ai+3, . . . ,n  and  columns  l,2,...,i+l  are  zero. 
Columns  r+2 , r+3 , . . . ,n  are  equal  and  are  partitioned  from 
the  remainder  of  the  matrix. 

It  follows  from  (3*33)  that 


4 


. 
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(3.41)  |Q^1)-Q^,1)|  <  |Q1||Q2|  ...  I  Qr_2  I  I E!  I 

+  I  Q  I  I  Q  2  I  •••  I  ®jr-3  11^2 

+  ...  +  |  Q-^  |  |  Q2  |  |  I 
+  1^1 1  lEr_2 I  +  I Er-1 1 


By  substituting  (3.40)  in  (3.41)  for  i=l f 2  , . . .  ,r-2 , 
obtain 


we 


(3-42) 
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Let  the  matrix  on  the  right-hand  side  of  (3.42)  be  denoted 
by  H  =  (h. . ) .  Then 

-L  J 


(3.43) 


h.  . 
ij 


(j-i+1)  2j  1  2  e, 
•  (r-i+1)  2r_i_1  e, 
0  ,  otherwise. 


From  (3.42),  we  have 


i  <  j  £  r 

j  >  r  and  i  <  r-1 


(3.44) 


< 


/0r-2  Is 
r ( 2  --)  + 


2r  ^  r(n-r+l) 


(n-r)r2r  ^ 


and  therefore,  can  state 

Theorem  3 . 4  Let  Fj  be  the  matrix  of  round-off 
errors  in  the  calculation  of  Q  in  (3.9).  Then 


(3. 45)  IIfJL  <  2r  2  r(n-r+l)  e. 

This  bound  is  at  least  as  small  as  that  given  by  (3.38) 
since 


2(n-r+l)  <  2n-r 


57 


and 


r  <  n(r-l) 


for  r  >  2. 

3.2.3  A  Bound  on  ll^loo  we  denote  the  exact  and 


computed  values  of  by  A^,1^  and  A^1^ 

then  Wilkinson  has  shown  in  [21]  that  A^1^ 
computational  equation 


,  respectively 
satisfies  the 


0.^6) 


j 


where 


. 


(3. 47) 
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where  the  left  partition  of  the  matrix  in  (3.47)  is  square 
and  of  order  m  and  where  g  is  the  element  of  maximum 
modulus  at  any  stage  in  the  reduction  of  A  to  A^\ 

The  only  round-off  errors  incurred  in  the  computation 
of  A  are  in  the  calculation  of  the  elements  d^ ,  i=l,2,. . . ,r, 
since  all  other  elements  of  A  may  be  set  to  zero.  If  A 
satisfies  the  computational  equation 


(3.48) 


A  =  A  +  5 


then  it  follows  from  (3.47)  that 


(3.49)  |E2|  <  ( 2 . 01 ) ge 


0  0  0  0 

0  1  0  .  .  .  0 

o 

• 

• 

• 

C\l 

o 

o 

0 

i — 1 

•  *  i 

u 

• 

• 

•  •  o 

•  •  o 

■  •  o 

0 

_ 

__ 

where 


E 


3 


is  a  diagonal  matrix  of  order  r. 


i  - 
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Let  e^,  i=l,2,...,r,  be  the  diagonal  elements  of  E^, 


i .  e . 


(3.50) 


e.  =  (i-1) ( 2 . 01 ) ge  . 


Therefore,  if  we  denote  the  exact  and  computed  values  of  the 
nonzero  terms  of  A  by  (d^)^  and  (d^)^,  respectively, 
then 


(3.5D 


(d±)E  -  (d± )  c  I  s  ,  1=1,2,. ...r 


From  the  computational  equation 


(3.52) 


d .  =  d .  +  e  .  , 

i  i  i  5 


we  obtain 


(3.53) 


1  _  1 


d .  d . +e . 

i  li 


+  e 


(1) 


where  it  is  assumed  that  d^+e^  /  0,  and  e|^  represents 
the  round-off  error  due  to  a  floating-point  division.  Hence, 
by  Theorem  2.1, 


(D  = 

i  d^+e 


(3.54) 
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We  may  rewrite  (3.53)  as 


(3.55) 


J_  =  JL  ei  .  -CD 

d.  "  d.  d. (d.+e.  )  i 

i  i  iii 


Therefore,  denoting  the  exact  and  computed  values  of  — 


by  (^p-)g  and  (^-)p,  respectively,  we  have  that 


d.  'C 

l 


(J_)  _  (_L) 

vd.;E  vd.;C 

l  l 


e . 

£  1 

d.+e.  -  d. (d.+e. ) 
i  l  l  i  i 


d .  e  -  e  . 
i _ i 

d. (d. +e .  ) 

l  i  l 


Thus,  if 


(3.56) 


(A-)'  e  (A-)’  +  P2  , 


then 


(3.57) |F2 
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The  result  of  this  section  is  stated  as 

Theorem  3 . 6  Let  F ^  be  the  matrix  of  round-off 
errors  in  the  oaloulation  of  ( )  '  3  defined  in  (Z.5).  Then 


(3.58) 


max 
1  <i  <r 


d  .  e  -  e  . 

V  X, 

d  . ( d  ,+e  . ) 

t  X,  X 

- 

3.2.4 
values  of 


A  Bound  on  II^Hoo  Let  the  exact  and  computed 
P,  defined  in  (3.8),  be  Pp  and  Pn.  Then 


(3.59) 


Since  the  matrices  FL  ,  i=l,2,...,r,  are  permutation 
matrices,  there  is  no  round-off  error  introduced  due  to  a 
matrix  multiplication  involving  any  of  the  FL  (see  Section 
2.3). 

Let  the  matrices  S^  be  defined  recursively  as  follows 


(3.60) 


P  R 
r 


r 


S.  =  S .  nP  .  ,-,R  .  , , 

l  l-l  r-i+1  r-i+15 


i=2 , 3, . . . ,r . 


Then 


' 
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(3.61) 


S  =  P  . 
r 


An  upper  bound  for  the  round-off  error  incurred  in 
the  matrix  multiplication  si_ipr_i+iRr_i+i  is  given  by 
(2.47).  If 


(3.62)  S.  ,P  .  ,R  .  n  =  S.  ,  P  .  ,,R  .  +  E  .  ,,  . 

l-l  r-i+1  r-i+1  l-l  r-i+1  r-i+1  r-i+1  5 


then 


(3.63)  II E  . 


<  S.  ,  P  .  ,  ,  R  . 
r-i+1"  °°  11  l  —  l*  oo  I*  r-i  +  1  r-i+11'00 


e ,  i=2 ,3, . . • ,r . 


From  (3.60)  and  (3.62)  we  obtain  the  following  set  of 
computational  equations  for  the  calculation  of  Pc : 


(3.64) 


S,  =  P  R 
1  r  r 


S0  =  P  R  P  nR  +  E  , 
2  r  r  r-1  r-1  r-1 


PRP  nR  -.P  0R  0  +  E  P  0R  o  +  E  0 
r  r  r-1  r-1  r-2  r-2  r-1  r-2  r-2  r-2 


Sr  E  PrRrPr-lRr-l* ’ * P1R1  +  Er-lPr-2Rr-2Pr- 3Rr- 3  *  *  * P1R1 

+  E  oP  -)R  oP  I,  R  I,  •  •  •  Pn  Rn  +  ...  +  E0P0R0P-|  R, 
r-2  r- 3  r-3  r-4  r-4  11  32211 


v. 


+  ^2P2Rl  +  Rl  * 


Hence,  from  the  final  computational  equation  in  (3.64), 
we  obtain 
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(3.65) 


|  Fj  <  ||  E  .P  0R  0P  0R  0  .  .  .  Pn  R,  „ 

'I  3"  oo  H  r-1  r-2  r-2  r-3  r-3  1  1  00 


+  ||  E  0P  0R  0P  i.  R  i.  •  •  •  P-,  R, 
11  r-2  r-3  r-3  r- 4  r- 4  11 


00 


+  ...  +  )l  E  2  P  2  R  2  ^  1  ^  1  ^  oc 
+  I|e2p1r1||oo  +  II E-JI  m  . 


Using  the  property 


(3.66) 


l|AB|L  *  ||A||J|B|| 


and  the  fact  that 


(3.67) 


I!  Rill  00 


=  1 


1-1,2,. . . ,r  1, 


the  inequality  (3.63)  becomes 


r-i  +  1 


II  Si_ i II  coll  Pr-i  +  lH  °°e 


< 


II  Prll  coll  Pr- 


r-i  +  2 


Pr-i+l 


(3.68) 
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Using  (3.66),  (3.67)  and  (3.68),  (3.65)  becomes 
(3-69)  ||P3IL  s  (r-l)||Pr||JPr_1IL---||P1ILe  . 

Since  all  of  the  multipliers  used  in  the  construction  of 
the  matrices  ,  1=1, 2,..., r,  have  the  property  that  their 
moduli  are  less  than  or  equal  to  one,  we  obtain 

(3.70)  flPjL  <  2  . 

The  following  theorem  is  a  result  of  (3*69)  and  (3.70). 

Theorem  3.6  Let  F  be  the  matrix  of  round-off  errors 
in  the  calculation  of  P  according  to  equation  (3.8).  Then 

(3.71)  ||P3IL  *  (r-l)2rz  . 

3.2.5  A  Bound  for  the  Computed  Reflexive  Generalized 
Inverse  A  norm  bound  for  the  round-off  error  incurred  in 
computing  a  Reflexive  Generalized  Inverse  can  be  obtained 
from  (3.16).  It  follows  from  (2.47)  that 

(3.72)  II^IL  5  e||Q|IJ(A')'IL 


and 
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(3.73)  |P5IL  ^MqIJ(a“)1J|pIL  • 

Substitution  of  (3.45),  (3.71),  (3.72)  and  (3.73)  into 
(3.16)  gives  the  following  norm  bound: 


<3.74)  ||  Ge  -  Gclco 


r (n-r+1 ) 2 


+  IIqIIJf 


r-2 


00 11  2  00 


+  II  p2IIJe 


P|l  +  (r-l)2  e 


00 


x( II  (a  )  ML  +  ||f2IU 


00 


llQlL  +  r (n-r+1 ) 2r  2e 


+  2e||Q||J|(A“)MIJ|P|L  > 


where  ||  F^ll  00  b°unded  in  (3.58).  The  product  II  P^||  J|  P3II  ^ 

2 

is  neglected  as  it  is  of  the  order  e  .  Our  computer  program 
for  the  computation  of  the  Reflexive  Generalized  Inverse 
includes  the  evaluation  of  the  right-hand  side  of  (3.74) 

(see  Chapter  V) . 


i  ■  ■ 


CHAPTER  IV 


A  NORM  BOUND  ON  THE  ERROR 
IN  COMPUTING  THE  PSEUDOINVERSE 

4 . 1  Computation  of  the  Pseudoinverse 

We  will  compute  the  Pseudoinverse  A+  of  an  arbitrary 

complex  matrix  A  of  dimension  m-by-n  according  to  the 

algorithm  given  by  Willner  [23].  Throughout  this  chapter, 

it  is  assumed  that  the  rank  of  A  is  r  and  that 

A.,  j=l,2,...,n,  is  the  j-th  column  of  A.  Also,  let 
J 

e^,  i=l , 2 , . . . ,m,  be  the  i-th  unit  vector  of  the  identity 
matrix  of  order  m.  The  proof  of  Winner’s  algorithm 
requires  the  following  lemma,  which  is  due  to  Greville  [7]. 

Lemma  4 . 1  If  a  matrix  A  of  dimension  m-by-n 
and  rank  r  can  be  expressed  as  a  product 

A  =  PB  , 

where  P  and  B  are  matrices  of  dimension  m-by-r  and 
r-by-n>  respectively ,  and  both  are  of  rank  r s  then 

(4.1)  A+  =  B*(BB*)~2 . 

Proof :  To  prove  the  above,  it  suffices  to  show 

that  equations  (1.1),  (1.2),  (1.3)  and  (1.4)  are  satisfied. 


■ 
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AA+A  =  PBB*(BB*)  1(P*P)“1P*PB 


=  PB 


=  A  . 


A+AA+  =  B*(BB*)“1(P*P)"1P*PBB*(BB*)"1(P*P)_1P* 
=  B*(BB*)"1(P*P)“1P* 


=  A+  . 


( A+A) *  =  [B*(BB*)_1(P*P)_1P*PB]* 


[B* (BB* )”1B] * 
B*[ (BB* )”1]*B 


B*[(BB*)*]“1B 

B*(BB*)_1B 

B*(BB*)“1(P*P)_1P*PB 


=  A  A  . 


( AA+ ) *  =  [PBB*(BB*)“1(P*P)"1P*]* 
=  [P(P*P)_1P*]* 

=  P[(P*P)-1]*P* 

=  P[ (P^P)*]’1?* 

=  P(P*P)-1P* 

=  PBB*(BB*)"1(P*P)"1P* 


=  AA 


Q .  E .  D . 
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Theorem  4.1  Let  A 
of  dimension  m-by-n  and  rank 
Eermite  normal  form  of  A}  and 
matrix  of  order  m  satisfying 


be  an  arbitrary  complex  matrix 

be  the 

let  Q~2  be  the  nonsingular 


r. 


Let  H  = 


B 

0 


(4.2) 


H  =  Q~2A  . 


If  Q  is  partitioned  as  Q  =  [P|i?],  where  P  is  of 
dimension  m-by-rs  then 


(4.3) 


A+  =  B* (P*AB* )~2  P*  . 


Proof :  Since  A  is  of  rank  r,  it  follows  that 

B  is  of  dimension  r-by-n.  The  existence  of  the  nonsingular 
matrix  Q-'1',  satisfying  (4.2),  follows  from  the  fact  that  A 
can  be  reduced  to  H  by  elementary  row  operations.  From  (4.2), 


A  =  QH 


=  [P|R] 


B 

0 


=  PB  . 


As  A,  P  and  B  are  all  of  rank  r,  we  can  use  (4.1)  to 
compute  A+.  Therefore, 
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A+  =  B*(BB*)"1(P*P)"1P* 

=  B^CP^PBB*)"1?* 

=  B* (P*AB* )_1P*  . 

Q  .  E  .  D . 

Winner’s  Algorithm  Willner  summarizes  the  compu¬ 
tation  of  A+  using  (4.3)  in  the  following  five  steps: 


(i) 

Compute  H  by  using  Gauss-Jordan 

transformations . 

(ii) 

From  H  determine  P  as  follows 

:  the 

i-th 

column  of  P,  P  ,  i=l,2,...,r,  is 

equal 

to  A . 

J 

if  Hj  =  e±s  j  =  1 , 2  ,  .  .  .  ,  n . 

(iii ) 

Calculate  P*AB*. 

(iv) 

Invert  P*AB*. 

(v) 

Calculate  A+  using  (4.3). 

We  now  proceed  to  place  an  upper  bound  on  the  round-off 
error  incurred  in  computing  A+  using  this  algorithm. 

4 . 2  An  Error  Analysis 

4.2.1  A  Bound  for  the  Computed  Hermite  Normal  Form 
The  algorithm  due  to  Marcus  and  Mine  (see  Section  2.1)  for 
the  computation  of  the  Hermite  normal  form  of  A  employs 
Type  II  elementary  row  operations.  As  the  elements  both 
above  and  below  the  pivotal  element  are  reduced  to  zero. 
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these  elementary  row  operations  are  equivalent  to  a  Gauss- 
Jordan  transformation,  i.e.,  pre-multiplication  of  A,  or 
a  reduced  form  of  A,  by  a  matrix  as  defined  by  (2.23). 

As  A  is  of  rank  r,  r  Jordan  transformations  are  required 
to  reduce  A  to  its  Hermite  normal  form.  As  in  Section  2.1, 
we  will  let  n^,  i=l,2,...,r,  denote  the  column  numbers  of 
the  r  columns  of  A  which  contain  pivotal  elements  for 
the  Jordan  transformations.  These  are  not  necessarily  the 
first  r  columns  of  A.  If  we  denote  the  matrix  A  after 
the  i-th  Jordan  transformation  by  A^1^  =  (afj^),  then  the 
multipliers  used  to  construct  J^  will  be  defined  by 


(4.4) 


a 


a 


(i-1) 
j  >n± 

(i-1)  5 

i,n± 


i  1 , 2  , .  .  .  ,  r  , 
j— lj2,...,m, 


rather  than  by  (2.24).  It  is  assumed  that  A ^  ^  =  A. 

As  all  of  the  pivotal  elements  of  the  Hermite  normal 

/  •  \ 

form  are  equal  to  one,  the  i-th  row  of  A ' 1  must  be  divided 

(  #  ^ 

by  the  pivotal  element  a.  .  This  may  be  accomplished  by 

l ,  n^ 

pre-multiplication  of  A ^ 1 ^  by  a  diagonal  matrix  =  (djk). 


where 


- 
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(4.5)  D±  = 


a. 


1 

nr 


ijni 


i 


i 


and 


d.  .  = 

li 


a ; 


1 

(IT 


If  we  let 


(H.6) 


Ja)  = 


DiJi 


i 


then 


(4.7)  J 


(1) 


1  0  ...  0  m, . 

li 


0  1  ...  0  m^ 


0  0 


0 


a . 


1 

nr 


i»n. 


0  0  .  .  .  0 


m  . 

mi 


0  ...  0 


0  ...  0 


0  0  ...  1  m.  ,  .  0  ...  0 

i—i,i 


0 


0 


0  0  ...  0  m. . i  .  1  ...  0 

i+l,i 


0  ...  1 


Therefore,  if  H  denotes  the  Hermite  Normal  form  of 


then 


(4.8) 


H  =  D  J  D  nJ  , 
r  r  r-1  r-1 


ViA 


T(i)Td) 

u  „  -i  •  •  • 

r  r-1 


The  intermediate  matrices  are  given  by 


(4.9) 


A 


(i)  _ 


T(1)T(1) 

l  l-l 


J^A 


j(Da(1  1) 

l 


i-1,2,.. 


. 


7^ 


and 

(4.10)  A^r)  =  H  . 

Complete  pivoting  may  not  be  used  since  H  must  be 
obtained  from  A  by  elementary  row  operations  only. 

Partial  pivoting,  however,  may  be  used  to  insure  that  those 
elements  of  below  the  main  diagonal  are  less  than 

or  equal  to  one  in  magnitude. 

/  .  \ 

It  is  possible  to  multiply  the  pivotal  row  of  A'  , 

i=0 , 1 , . . . ,r-l ,  by  a  power  of  2  (in  a  binary  machine)  to 

insure  that  all  elements  of  will  be  less  than  or 

l+l 

equal  to  one  in  magnitude.  We  let  P  ,  i=0 , 1 , . . . ,r-l ,  be 
equal  to  a  unit  matrix  except  for  its  (i+l)-th  diagonal 
element,  which  is  some  non-negative  power  of  2,  i.e. 

1 

1 

1 

( i ) 

pi+l,i+'l 

1 


(4.11)  P±  = 


1 


.  • 
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where 


(4.12) 


( i ) 

^1+1 ,1+1 


is  chosen  so  that  the  modulus  of  the  product  of  2^i  and 

/  •  \  /  •  \ 

the  pivotal  element  of  A^1  ,  a.^n'  ,  is  greater  than  or 

,ni  +  l 

equal  to  |a.^  |,  j=l,2,...,m,  j^i+1.  Then 

J  >ni  +  1 


(4.13) 


H  =  J(1)P  ,J(1)P  ? 

r  r-1  r-1  r-2 


J  ^  ^  P  A 
J1  rCT 


9 


where  all  elements  of  each  matrix  1=1, 2,..., r,  are 

less  than  or  equal  to  one  in  magnitude. 

Since  llP-JI^  Is  unbounded,  this  does  not,  however, 
lower  the  upper  bound  for  the  round-off  error  in  the  calcula¬ 
tion  of  H.  Hence,  we  will  use  (4.8)  for  the  computation  of 
H. 

If  we  let  ,  i=l,2,...,r,  be  the  matrix  of  round-off 
errors  incurred  in  multiplying  A^  ^  by  jf1),  then  the 
computed  value  of  H,  Hc,  is  defined  by  the  following  set 
of  computational  equations : 


•  1  s  i+i,i+i-q 


.  S*!*!  ,  T»  lo  £:^neTi9l9  IIb  9*i9riw 

.  i-  OJ  .  0  .  UC  : '0  -J  *io 


'• 


76 


(4.14) 


A(1)  E  J^A  +  E1 


A(2)  E  J^J^A  +  J^1^!  +  e2 


A(r)  E  jCDjCl)  .  .  .  J(1)A 
r  r-1  1 


,  t(dt(i)  T(i)F 

r  r-1  2  1 


+  J(1)J(1]  ...  J^Ep  +  .  .  . 
r  r-1  3  2 


+  J(1)E  n  +  E 
r  r-1  r 


Since  (4.8)  defines  the  exact  value,  Hg,  of  H,  it 
from  the  last  equation  in  (4.14)  that 


(4.15) 


|He-hc|.  =  •••  J2";ei 


(i). 


+  ...  J^E  +  . 

r  r-1  3  2 


+  J(1)E  ,  +  E 

r  r- 1  r" 00 


*  Il41)l  JJr-lL  •••  IUp1} 


+  lJr1)LIJr-iL  •••  ll^15 


(1) 


+  •••  +  lJr  LlVlL  + 


follows 


JEJ. 

IIJ|E2|| 

«ErIL 


.<  vtfLftv  r- .©*  -  <  .  e  13 

3  ;rftf  0  .  n  :  rol^fijjpa  l  urfl  rtioT: 
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By  (2.47), 


(4.16) 


II E,  || 


<  e  J 


(1) 


II  A 

00  M 


(i-1) 


00 


i=l ,  2 ,  .  .  .  ,r 


Therefore,  using  (4.9), 


(4.17)  ||  E±| 


00 


<  e 


Ja) 

i 


1  =  1,2 


Substituting  (4.17)  in  (4.15),  we  obtain 


(4.18) 


Thus  we  may  state 

Theorem  4.2  A  norm  bound  on  the  round-off  error 
incurred  in  computing  the  Hermite  normal  form  of  A  is 


where  r  is  the  rank  of  A  and  the  matrices 
defined  by  (4.7). 


Jn) 


are 


4,2.2  A  Bound  for  the  Product  P*AB*  The  determination 
of  the  matrix  P  by  Winner’s  algorithm  involves  no  arithmetic 
calculations.  Therefore  no  round-off  errors  are  introduced. 


■ 


Ol  .  .►vXrWf  t  t£  ,  .  i;.  a1  V  a  X±t^j5i'T  fJ.ij  lo 


tb90jjbort^n  sr  s  r  s  :o  *10-  ■rjc-  or  vu  >ii  J  .  snoldBluoXso 


78 


and  the  elements  of  P  are  exact  (to  the  precision  of  the 
computer ) . 

As  the  initial  step  in  bounding  the  round-off  error 
incurred  in  computing  the  matrix  product  P*AB*,  we  consider 
the  computation  of  P*A.  If  is  the  matrix  of  round-off 

errors  for  the  multiplication  of  P*  by  A,  then 


(4.19) 


P*A  =  P*A  +  F  , 


and  by  (2.47), 


(4.20) 


e 


00 


The  elements  of  B*  are  not  exact,  as  are  the  elements 
of  A  and  P*,  because  B  is  a  submatrix  of  H.  It  follows 
from  (4.14)  that 


HS-H8 


(Djd) 

r  r-1 


E 


1 


+  J 


ayi) 

r  r-1 


J(1)E  , 

r  r-1 


E  )* 
r 


or 


. 
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(4.21) 

|H*-H*|  =  |EfJ^1)*J^1)*  ...  J^1}* 

+  E-J*1**^1’*  ...  j<1>*  +  ... 

+  +  E*l  • 

Therefore , 


(4.22) 

C\ J 

+ 

M 

VI 

*o 

where 


(4.23) 

F2  =  Ie-j^'j^*  ... 

+  E*J i1 ^  *J  h1 ^  *  ...  J^*  +  ... 

2  p  4  r 

+  E*  J  ^ 1  ^  *  +  E*  1  . 

r-1  r  r1 

Thus 


(4.24) 

lH*IL  s  »h|IL  +  II  p2IL  • 

We  will  denote  the  exact  and  computed  values  of  B  by 
Bg  and  Bc,  respectively.  Since  B*  is  a  submatrix  of 


H§»  and  ||  H||  „  =  ||B*||ro, 


it  follows  that 


*****  Bwoixcrt  it  ,J§aj  -  „|*H|  bns  ,.»H 
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(4.25) 


Halloo  *  IB|IL  +  I|p2IL  . 


The  corresponding  computational  equation  is 


(4.26) 


B*  =  B*  +  PQ  , 


where 


(l)*n  n  T ( 1 ) * 


(4.27)  ||p2L  s  lE*ILII J2  ;  II Jj 


00 


•  ••  II 


* 


+  |E?IJji1)*|Jji1)*|0„  ...  |J<1)-| 


+  •  •  •  +  II  E*_ill  oo I  Jr 


(l)* 


I  E?H  00 


=  lEllllJ21)lllJ31)ll  •••  llJr±;|ll 


(1) 


llE2llillJS1)llillJ41)llx  •••  HJrx;H 


(1) 


+  ...  +  ||Er_i||  x||  J^1)||  1  +  II  E|| 


s  rell  )|l  ].ll  Jr-lll  1  •••  llJn(1) 


1  "  1  * 


The  computational  equation  for  the  product  P*AB*, 


using 


(4.19)  and  (4.26),  can  now  be  written  as 
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(4.28)  P*AB*  E  (P*A+F1) (B*+F2)  +  F^ 


or 

(4.29)  P*AB*  E  P*AB*  +  P*AF2  +  F1B*  +  F^  +  F3  , 

where  F^  Is  the  matrix  of  round-off  errors  incurred  in  the 

multiplication  of  P*A  by  B*.  The  term  F^F2  and  the 

round-off  error  incurred  in  forming  the  products  P*AF2  and 

2 

F-^B*  will  be  neglected  since  they  are  of  the  order  e  . 

By  (2.47), 

(4-30)  ||F3IL  s  e||P*A||J|B*||0o  . 

The  following  theorem  is  a  result  of  (4.29)  and  the  above 
discussion . 


Theorem  4.3  If 


(4.31) 


V'AB*  E  P*AB*  +  F  , 


then 


(4.32) 
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where  the  hounds  on  ||^H<»*  ||  oo  an ^  ll^^lloo  ave  9^ven  by 
(4.20)3  (4.27)  and  (4.30)3  respectively. 

4.2.3  A  Bound  in  Computing  (P*AB*)~1  The  fourth  step 
in  Winner’s  algorithm  for  the  computation  of  A+  is  the 
inversion  of  the  square  matrix  P*AB*  of  order  r.  For 
the  computed  inverse,  X”1,  °f  a  nonsingular  matrix  X 
of  order  n,  Wilkinson  [21]  has  shown  that 

(4.33)  XX"1  -  I  =  K  , 

where 

(4.34)  II  K||  *  ||  X"1!^  g(2.005  n2+n3)£  , 

2 

with  terms  of  order  e  ignored.  g  denotes  the  element 
of  maximum  modulus  at  any  stage  in  the  reduction  of  X  to 
X"1.  Therefore, 

(4.35)  (P*AB*+F) (P*AB*+F)”1  -  I  =  K  , 

where 


(4.36) 


|j  K|| ^  <  II  (P*AB*+F)c1||ro  g(2.005r2+r3)e  . 


- 
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Hence 


(4.37) 


(P*AB*+F)C1 


(P*AB*+F)"1(I+K)  , 


where  (P*AB*+F)g^  denotes  the  exact  inverse  of  (P*AB*+F). 

Wang  [19]  has  shown  that  the  inverse  of  a  modified 
matrix  Y+E  can  be  obtained  from 


(4.38) 


(Y+E)  1  =  ( I+Y  1E )  1Y_1  . 


Using  (4.38)  in  (4.37),  we  obtain 

(4.39)  (P^AB^+F)"1  =  [I+(P*AB*)g1F]"1(P*AB*)g1(I+K) 

=  ( P*AB* ) +  [I+(P*AB*)”1F]”1(P*AB*)g1(I+K) 

-  ( P*AB* )”1 
E 

=  (P*AB*)“1  +  [I+(P*AB*)"1F]_1(P*AB*)g1 
x [K-F(P*AB*)“1] 

=  ( P*AB* ) +  (P*AB*+F)"1[K-F(P*AB*)”1] . 

H  Hi  H 


From  (4.39),  we  can  write  the  computational  equation 
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(4. 40) 

where 

(4.41) 

Prom  (4.39), 

(4.42) 

Substituting 

(4.43) 

where 


(P*AB*)C1  S  (P*AB*)"1  +  F4  , 

l|F4L  £  ||  CP*AB*+P)-1|[oo||  K-F(P*AB*)~1||  ^  . 
it  also  follows  that 
||  (P*AB*  +  F)‘1||oo  <  ||(P*AB*)-1||ro 

+  ||  (P*AB*  +  F)"1||  J|K-F(P*AB*)g1||oo. 
(4.36)  in  (4.42),  we  obtain 
||(P*AB*  +  F)-1||oo  <  ||(P*AB*)-1||oo 

+  ||  (P*AB*  +  F)~1||  ^  TcH  (P*AB*+F)”1||ct 

+  |F|IJI(P*AB*)'1|I00J, 

c  =  g(2 . 005r2+r^)e  . 


(4.44) 


4 
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Assuming  that  the  elements  of  |P|  are  small  relative  to 
those  of  |P*AB*| ,  and  that  the  matrix  P*AB*  is  well- 
conditioned,  the  exact  and  computed  values  of  ||  (P*AB*+F)”'*"||  a 
are  approximately  equal.  In  this  case,  (4.43)  becomes 


(4.45) 


—  c ||  (P*AB*  +  F)”1||^  + 


1-ll  Fl  ool  (p*ab*)e1IL 


X  ||  (P*AB*+F)  1||oo  -  ||(P*AB*)e1||co  <  0  . 


The  discriminant  of  this 
||  (P*AB*+F)_1||  is 


quadratic  expression  in  the 


variable 


(4.46)  d  =  (1-11  F||  J|  (P*AB*)e1||co)2  -  4c||  (P*AB*  )E1|| 


which  we  can  assume  to  be  positive  since  c  is  of  the  order 
e.  It  follows  that  the  roots  are  either  both  positive  or 
both  negative,  since  the  value  of  the  quadratic  expression 
is  negative  for  ||  (P*AB*+F)_1||  ^  =0.  To  place  an  upper 
bound  on  ||  (P*AB*+F)_1||  ,  using  (4.45),  we  need  consider 

only  the  case  of  two  positive  roots.  For 


(4.47) 


1  -  II  FIJI  (P*AB*)”1|L  >  0  , 


both  roots  will  be  positive,  and  the  smaller  of  the  two  roots 
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will  equal 


(4.48) 


-1  +  II  F|  Jl  (P*AB*)”1||co  +  d1/2 


-2c 


Then  (4.45)  implies  that 


-In  .1/2 


(4.49)  II  (P*AB*+F)  X||  < 

■ 1  00 


1  -  ||  F||J|  (P*AB*  )-x||  -  d 


2c 


and  from  (4.36)  and  (4.4l),  we  have  that 


(it. 50)  ||F4||00  s  ||  (P*AB*+F)_1|| 


+  ||  F||  Jl  ( P*AB*  ) 


c||(P*AB*  +  F)  1|| 


CO 


-ll 


The  result  of  this  section  is  summarized  in  the  following 


theorem. 


Theorem  4.4  Given  that 


P*AB*  =  P*AB*  +  F  , 


where  ||  F||  „ 


is  bounded  by  (4.32),  then 
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(ptAB*)'1  E  (P+AB*)'1  +  F4  3 

where  is  bounded  by  (4.50). 

4.2.4  A  Bound  for  the  Computed  Pseudoinverse  The 

final  step  in  Winner's  algorithm  for  the  computation  of 

+  I  _  i 

A  is  the  calculation  of  the  matrix  product  B*(P*AB*)  P*. 

Considering  first  the  multiplication  of  B*  by  (P*AB*)”1, 
it  follows  from  (4.26)  and  (4.40)  that 

(4.51)  B*(P*AB*)-1  e  (B*+F2)[(P*AB*)“1  +  F^]  +  F^ 
or 

(4.52)  B*(P*AB*)-1  5  B*(P*AB*)-1  +  B*Fi)  +  F2(P*AB*)-1 

+  f2f4  +  f5  , 

where  F._  is  the  matrix  of  round-off  errors  in  the  multi- 
5 

plication  of  B*  by  (P*AB*)-1.  Since  the  product  F2F4 

and  the  round-off  errors  incurred  in  forming  the  products 

-1  2 

B*F^  and  F2(P*AB*)  are  of  order  e  ,  they  are 
neglected.  Therefore, 

B*(P*AB*)_1  =  B*(P*AB*)-1  +  B*F4  +  F2(P*AB*)_1  +  F  , 


(4.53) 


■ 
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where ,  by  (2.47), 


(“•54)  ||F5IL  s  e||B*||  J  (P*AB*)-1||0o  . 

Post-multiplication  of  (4.53)  by  P*  produces 
(4.55)  B*(P*AB*)"1P*  E  B*(P*AB*)-1P* 

+  B*F^P*  +  F2(P*AB*)_1P* 

+  p  p#  +  p 

5  6  5 

where  Fg  is  the  matrix  of  round-off  errors  incurred  in 
the  multiplication  of  B*(P*AB*)-1  and  P*.  Hence 

(“•56)  llPgIL  £  e||B*|IJ|(P*AB*)-1||Jp*||oo  . 

Since  A+  =  B*(P*AB*)  it  follows  from  (4.55)  that 

(4.57)  II  AE-AclL  s  |b*IJ|p4|Jp*I<» 

+  IIf2II  J  (p*AB*)-1ll  Jp*IL 
+  I|f5IIJ|p*IL  +  IlFglloo  . 
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where 


A+  and  A+ 
respectively . 


denote  the  exact  and  computed  values  of 
By  (4.54)  and  (4.56), 


(4.58) 


The  result  of  substituting  (4.50)  and  (4.58)  into  (4.57) 
is  stated  in  the  following  theorem. 

Theorem  4 . 5  If  the  Pseudoinverse  A+  of  an 
arbitrary  matrix  A  is  computed  using  Willner's  algorithm 
(see  Section  4.1) ,  then  the  round-off  error  in  the  computed 
value  of  A+  is  bounded  by 


(4.59)  II £  II p*l 


IMI JF 


ool'  ^  00 


-1 


+  II  ( P  *AB*  )  'i|l00(|l^2ll00^2e||B^|oo; 


where  F ^  and  F are  defined  by  (4.26)  and  (4.40) , 
respectively ,  and  they  are  bounded  by  (4.27)  and  (4.50), 
respective ly . 

||  B * ||  ^ ,  ||  P * ||  oo  and  ||  ( P*AB#  )-1||  ^  must  be  computed  in 
the  process  of  calculating  A+ .  Our  computer  program  for 
the  computation  of  the  Pseudoinverse  includes  the  evalua¬ 
tion  of  the  right-hand  side  of  (4.59)  (see  Chapter  V). 


CHAPTER  V 


NUMERICAL  RESULTS 

The  computation  of  a  Reflexive  Generalized  Inverse 
and  the  Pseuodinverse  was  performed  using  the  IBM  System 
360  FORTRAN  IV  source  language  and  the  IBM  System  360/67 
computer  system  at  the  University  of  Alberta.  Several 
test  matrices,  for  which  the  exact  Reflexive  Generalized 
Inverse  and  Pseudoinverse  are  known,  were  used.  The 
programs  evaluate  the  norm  bounds  for  the  round-off  error 
as  given  in  Chapters  III  and  IV,  and  a  comparison  is  made 
between  the  actual  and  predicted  errors. 

5 . 1  Test  Data 

5.1.1  Nonsingular  Matrices  The  test  data  for  both 
the  computation  of  a  Reflexive  Generalized  Inverse  and 
the  Pseudoinverse  are  based  on  a  set  of  nonsingular 
matrices  A^  of  order  n  with  known  inverses.  For  any 
positive  integer  k,  rectangular  matrices  of  dimension 
n-by-kn  with  known  Reflexive  Generalized  Inverses  and 
Pseudoinverses  can  be  constructed  by  catenating  k  such 
identical  nonsingular  matrices.  The  following  nonsingular 
matrices  denoted  A^ ,  A^ ,  A^  and  A^  are  due  to  Newman 
and  Todd  [12],  while  matrices  A^  and  A^  are  due  to 
Charmonman  and  Julius  [2], 


"  •  j| 
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The  matrix  is  defined  by 


(5.1) 


a .  . 

ij 


n+1 


1/2 


sin 


ij  tt 
n+1 


Since  is  symmetric  and  orthogonal,  we  have 


(5.2) 


A'1  =  A  . 


The  matrix  is  defined  by 


(5.3) 


a.  .  = « 
ij 


i/j  ,  i  -  j  ; 

j/i  j  1  >  J  3 


and  the  inverse  of  A^  is  a  symmetric,  triple-diagonal 


matrix  with 


(5.4) 


a.  . 

ij 


4j3 

4i2-l 


n 


2n-l 


-1(1+1) 

21+1 


0 


,  i=j  ,  i  <  n  ; 


•  • 


i=J=n; 


•  *1  n 

i-j I =i  ; 


i-j I >i  • 
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The  matrix  is  a  triple-diagonal  matrix  defined  by 


(5.5) 


a.  .  = 
ij 


-2  J  i=j  ; 


1  >  I i-j I =1  ; 


0  ,  | i-J  | >1  . 


Its  inverse  is  given  by 


(5.6) 


a .  . 
ij 


— 


-i  (n-,1  +1 ) 

n+1 


a .  . 

Ji 


i  -  J  j 


i  >  J 


The  matrix  is  defined  by 


(5.7) 


a  =  2  min  (i,j)  -  1  a 

J 


and  its  inverse  is  a  triple-diagonal  matrix  given  by 


(5.8) 


1.5 


i=j  =  l 


a.  .  =  i 
ij 


0.5  5 
1  , 
-0.5  , 


i=j=n  ; 
l<i= j  <n  ; 

I i-J l=i  ; 
I i-J |>i  . 


0 


) 
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The  following  definitions  are  used  to  define  the 

matrices  An  and  Ar  . 

5  6 

Definition  5 . 1  An  r-row  circulant  is  a  square 
matrix  of  order  n  in  which  the  i-th  row3  i=23Z3...3n3 
is  obtained  from  the  (i-l)-th  row  by  cyclically  shifting 
each  element  r  places  to  the  right. 

Definition  5.2  An  r- column  circulant  is  a  square 
matrix  of  order  n  in  which  the  i-th  column 3  i=23Z3...3n3 

is  obtained  from  the  (i-l)-th  column  by  cyclically  shift¬ 
ing  each  element  r  places  down. 

The  matrix  A._  is  the  r-row  circulant  with  first 

o 

row  (a,ah, . . . ,ahn  ,  where  a  and  h  are  arbitrary 

constants.  The  inverse  of  A,,  is  the  r-column  circulant 

5 

with  first  column  (b  ,0 ,0  , . . .  ,0 ,-hb) ' ,  where 


(5.9) 


b 


1 

a(l-hn) 


The  matrix  A^  is  the  r-row  circulant  with  first 
row  [a,a+h, . . . ,a+(n-l)h] .  The  inverse  of  A^  is  the 
r-column  circulant  with  first  column  (b-a ,b ,b , . . . ,b ,b+a) ’ , 
where 


2 


(5.10) 


b 


n2[2a+(n-l)h] 
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and 


(5.11) 


5.1.2  Rectangular  Matrices  with  Known  Reflexive 
Generalized  Inverses  Let  B  be  a  matrix  of  dimension 
m-by-n,  m<n,  with  a  nonsingular  matrix  of  order  m 

in  its  first  m  columns,  i.e. 


(5.12) 


5 


where  C  is  an  arbitrary  matrix  of  dimension  m-by-(n-m). 
It  can  easily  be  shown  that  a  Reflexive  Generalized 
Inverse  of  B  is 


(5.13) 


where  G  is  of  dimension  n-by-m.  Furthermore,  except 
for  round-off  error,  this  is  the  Reflexive  Generalized 
Inverse  that  is  calculated  by  the  elimination  method  given 
in  Section  3.1,  provided  that  each  pivotal  element  is 
chosen  only  from  the  first  m  columns  of  B.  This  latter 
condition  insures  that  rows  m+1 ,m+2 , . . . ,n  of  the  matrix 


t 
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Q  defined  by  (3.9)  are  unit  vectors.  Since  rows 
m+1 ,m+2 , . . . ,n  of  (A-)’  are  zero,  it  follows  that  rows 
m+1 ,m+2 , . . . ,n  of  the  product  Q(A-)’  must  be  zero.  This 
guarantees  that  rows  m+1 ,m+2 , . . . ,n  of  G  are  zero,  and 
therefore,  G  must  be  of  the  form  given  by  (5.13)  to 
satisfy  equations  (1.8)  and  (1.9). 

Complete  pivoting  is  used  in  the  Gaussian  forward 
elimination  for  the  reduction  of  the  matrix  B  in  (5.12). 
At  the  i-th  stage  in  the  reduction,  the  pivotal  element  is 
chosen  from  rows  i,i+l,...,m  and  columns  i,i+l,...,n. 
The  condition  that  the  pivot  be  chosen  only  from  columns 
i,i+l,...,m  will  be  satisfied  if  the  column  vectors  of 
C  are  equal  to  column  vectors  of  .  The  data  for  our 
program  has  been  constructed  by  setting  C  equal  to  A^ 
or  [A.  I  A. ]  ,  a  matrix  of  dimension  n-by-2n.  Therefore, 

B  is  of  dimension  n-by-2n  or  n-by-3n,  and  takes  the 
form 

(5.14)  CA± | A± ] 

or 

C  Ai  I  A±  I  Ai  ]  ’ 


(5.15) 


' 

■ 
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where  is  one  of  the  six  nonsingular  matrices  of  order 

n  defined  in  Subsection  5.1.1.  Then  the  Reflexive 
Generalized  Inverse  obtained  by  using  the  algorithm  of 
Section  3.1  is 


(5.16) 


a:1 

l 

0 


for  the  matrix  (5.14),  and 


(5.17) 


A 


-1 


0 


0 


for  the  matrix  (5.15).  The  dimensions  of  the  matrices 
(5.16)  and  (5.17)  are  2n-by-n  and  3n-by-n,  respectively 
Various  values  of  n  between  3  and  20  were  used  for 


computational  purposes. 


5.1.3  Rectangular  Matrices  with  Known  Pseudoinverses 
Let  B  be  a  matrix  of  dimension  m-by-km,  constructed 
by  placing  k  identical  nonsingular  matrices  A^ ,  of 
order  m,  side-by-side;  i.e. 


B  =  [A.  I  A.  I  A.  I  ...  I  A. ]  . 
i  1  i  1  i  1  1  i 


(5.18) 


« 
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Let  C  be  a  matrix  of  dimension  km-by-m,  constructed  by 
placing  k  matrices  cA”1  together,  where  c  is  a 
constant;  i.e. 


(5.19) 


Then 


(5.20)  BC  =  kcl  . 

Therefore , 

(5.21)  BCB  =  kcB 

and 

(5.22)  CBC  =  kcC  . 

For  C  to  be  the  Pseudoinverse  of  B,  equations  (1.1) 
and  (1.2)  must  be  satisfied,  i.e. 
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(5.23) 


BCB  =  B 


and 


(5.24) 


CBC  =  C  . 


From  (5.21)  and  (5.22),  it  follows  that 
(5.25)  c  =  i  . 


Since 


(5.26)  CB  =  i 


where  each  identity  matrix  is  of  order  m,  (5.20)  and 
(5.26)  verify  that  equations  (1.4)  and  (1.3)  are  satisfied. 
Therefore,  C  is  the  Pseudoinverse  of  B. 

The  data  for  our  program  have  been  constructed  by 
choosing  to  be  one  of  the  six  nonsingular  matrices 

defined  in  Subsection  5.1.1.  Only  values  of  2  and  3  have 
been  used  for  k,  with  n  assuming  various  values  between 


3  and  20. 


99 


5 . 2  Summary  of  the  Results 

The  Reflexive  Generalized  Inverse  given  by  (5.16) 
or  (5.17)  and  the  Pseudoinverse  given  by  (5.19)  and  (5.25) 
have  been  computed,  choosing  to  be  one  of  the  six 

nonsingular  matrices  of  Subsection  5.1.1.  The  matrix 
Aj-  was  generated  for  the  values  a=l,  h=1.5  and  r=l, 
and  a=l,  h=2  and  r=l  for  the  computation  of  the 
Pseudoinverse.  For  the  computation  of  a  Reflexive 


Generalized  Inverse, 

the 

values  a=l 

,  h=2  and 

OJ 

II 

or 

C\J 

ii 

i — 1 

ii 

cti 

and  r=3} 

and 

the 

values 

a=l ,  h=  3 

and 

r=4 

or  a=l , 

h=  3  and  r 

=  3 

were 

used  to 

generate 

V 

The 

second  of  each  pair  of  values  given  is  necessary  for  matrices 

Ar  of  certain  dimensions,  since  the  first  set  of  values 
5 

generates  a  matrix  of  rank  less  than  n.  The  matrix  Ag 
was  generated  for  the  values  a=l,  h=2  and  r=l,  and 
a=l,  h=3  and  r=l  for  the  computation  of  the  Pseudo¬ 
inverse.  For  computation  of  a  Reflexive  Generalized 


Inverse,  the  values 

i — 1 

II 

1 — 1 

II 

and 

r=2  or  a=l , 

ii 

i— 1 

and  r=3 ,  and  a=l. 

ii 

uo 

and 

r=4 

or  a=l ,  h=3 

and 

r=3  were  used.  Eight  different  nonsingular  matrices  were 
generated  for  the  construction  of  rectangular  matrices  of 
dimension  3-by-6,  3-by-9,  5-by-10,  5-by-15,  8-by-l6, 
8-by-24,  10-by-20 ,  10-by-30,  15-by-30,  15-by-45  and 
20-by-40. 


k  i  n  1  ■ 
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5.2.1  Results  for  the  Reflexive  Generalized  Inverse 
Table  5.1  gives  the  predicted  relative  error  and  the  actual 
relative  error  incurred  in  the  computation  of  the  Reflexive 
Generalized  Inverse  G,  given  by  (5.16)  or  (5.17),  of  a 
matrix  A  of  dimension  m-by-2m  or  m-by-3m.  The  pre¬ 
dicted  absolute  error,  ||Ge-Gc(|oo,  is  given  by  (3.7*0,  and 
the  predicted  relative  error  is  defined  to  be 


(5.27) 


The  actual  relative  error  can  be  computed  since  the  exact 

Reflexive  Generalized  Inverse  is  known. 

The  first  column  of  Table  5.1  lists  the  nonsingular 

matrices  of  Subsection  5.1.1  used  to  construct  the 

matrices  A,  which  take  either  the  form  of  (5.1*0  or 

(5.15).  The  three  parenthesized  numbers  after  A^  and 

Ag  are,  respectively,  the  values  of  a,  h  and  r  used 

to  construct  these  matrices.  The  second  column  of  the 

table  gives  the  dimension  of  A.  The  notation  aE±b 

+ 1^ 

denotes  the  value  axlO 
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TABLE  5.1 

ERROR  IN  THE  REFLEXIVE  GENERALIZED  INVERSE 


PREDICTED 

ACTUAL 

MATRIX 

DIMENSION 

RELATIVE  ERROR 

RELATIVE  ERROR 

Ai 

3-by-6 

0.7344E-4 

0.1641E-5 

A2 

3-by-6 

0.2532E-4 

0.7196E-6 

A3 

3-by-6 

0.4492E-4 

0.0 

A4 

3-by-6 

0.4475E-4 

0.5811E-7 

A5(l,2,2) 

3-by-6 

0.8338E-3 

0.0 

A5(l,3,^) 

3-by-6 

0 . 2739E-1 

0.1652E-6 

Ag  (1,1,2) 

3-by-6 

0.2941E-3 

0.4952E-6 

Ag(l,3,4) 

3-by-6 

0.1409E-1 

0.2682E-6 

A1 

3-by-9 

0.1075E-3 

0.1641E-5 

A2 

3-by-9 

0 . 3915E-4 

0.7196E-6 

A3 

3-by-9 

0 . 6689E-4 

0.0 

A4 

3-by-9 

0.6698E-4 

0 . 5811E-7 

A5(l,2,2) 

3-by-9 

0.1171E-2 

0.0 

A5(l,3,4) 

3-by-9 

0. 3860E-1 

0.1652E-6 

Ag(l,l,2) 

3-by-9 

0.4038E-3 

0.4952E-6 

Ag  ( 1 , 3, 4  ) 

3-by-9 

0.1923E-1 

0.2682E-6 

A1 

5-by-10 

0.7144E-3 

0.9780E-6 

A2 

5-by-10 

0.2189E-3 

0.4631E-6 

A3 

5-by-10 

0.1986E-3 

0. 33HE-6 

a4 

5-by-10 

0. 4468E-3 

0.1025E-5 

* 
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TABLE  5.1  (Continued) 

PREDICTED 


ACTUAL 


MATRIX 

DIMENSION 

RELATIVE  ERROR 

RELATIVE  ER 

(1,2,2) 

5-by-10 

0.5990 

0.0 

Ag(l,3,4) 

5-by-10 

0.37HE+3 

0.5297E-6 

Ag  ( 1 , 1 , 2 ) 

5-by-10 

0.8060E-2 

0.4572E-6 

A6(l,3,4) 

5-by-10 

0.5239 

0.4951E-6 

A1 

5-by-15 

0.1086E-2 

0 . 9780E-6 

A2 

5-by-15 

0 . 3464E-3 

0.4631E-6 

A3 

5-by-15 

0.3131E-3 

0. 33HE-6 

A4 

5-by-15 

0.6786E-3 

0.1025E-5 

A5( 1,2,2) 

5-by-15 

0.8273 

0.0 

A5  ( 1 , 3 , 4 ) 

5-by-15 

0 . 5301E+3 

0 . 5297E-6 

Ag(l,l,2) 

5-by-15 

0.1016E-1 

0.4572E-6 

A6(l,3,4) 

5-by-15 

0.6635 

0. 4951E-6 

A1 

8-by-l6 

0.1494E-1 

0.4615E-5 

A2 

8-by-l6 

0. 3796E-2 

0.1197E-5 

A3 

8-by-l6 

0.2662E-2 

O.3636E-6 

a4 

8-by-l6 

0 . 6602E-2 

0. 3008E-5 

Ag (1 , 2 , 3) 

8-by-l6 

0.4468E+4 

0.0 

Ag (1 j  3 , 3) 

8-by-l6 

0.1137E+8 

0.6392E-6 

Ag(l,l,3) 

8-by-l6 

0.1329 

0 . 7729E-6 

Ag ( 1, 3, 3) 

8-by-l6 

0 . 5950E+1 

0 . 37^9E-5 

A1 

8-by-24 

0.2484E-1 

0.4615E-5 

- 
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TABLE  5.1  (Continued) 


MATRIX 

DIMENSION 

PREDICTED 
RELATIVE  ERROR 

ACTUAL 

RELATIVE  ERROR 

A2 

8-by-24 

0.6180E-2 

0.1197E-5 

A3 

8-by-24 

0.4603E-2 

O.3636E-6 

8-by-24 

0.1037E-1 

0.3008E-5 

Ag(l,2,3) 

8-by-24 

0.6076E+4 

0.0 

A5(l,3,3) 

8-by-24 

0 . 1602E+8 

0 . 6392E-6 

Ag(l,l,3) 

8-by-24 

0.1608 

0 . 7729E-6 

Ag  ( 1, 3, 3) 

8-by-24 

0 . 7907E+1 

0.37^9E,-5 

A1 

10-by-20 

0.9892E-1 

0.8823E-5 

A2 

10-by-20 

0 . 2221E-1 

0.1199E-5 

A3 

10-by-20 

0 . 1648E-1 

0.3854E-6 

a4 

10-by-20 

0 . 3559E-1 

0 . 2895E-5 

Ag (1,2,3) 

10-by-20 

0.4699E+6 

0.0 

Ag (1 , 3, 3) 

10-by-20 

0.9169E+9 

0.7963E-6 

Ag(l,l,3) 

10-by-20 

0.3535 

0.7900E-6 

Ag (1,3,3) 

10-by-20 

0.2970E+2 

0 . 3735E-5 

Al 

10-by-30 

0.1713 

0.8823E-5 

A2 

10-by-30 

0.3665E-1 

0.1199E-5 

A3 

10-by-30 

0.2951E-1 

0 . 3854E-6 

A4 

10-by-30 

0 . 5765E-1 

0 . 2895E-5 

Ag (1,2,3) 

10-by-30 

0.6540E+6 

0.0 

Ag (l,3,3) 

10-by-30 

0.1310E+10 

0.7963E-6 
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TABLE  5 . 1 

(Continued ) 

MATRIX 

DIMENSION 

PREDICTED 
RELATIVE  ERROR 

ACTUAL 

RELATIVE  ERROR 

Ag(l,l,3) 

10-by-30 

0.4532 

0.7900E-6 

Ag  ( 1 , 3>  3) 

10-by-30 

0 . 3620E+2 

0 . 3735E-5 

A1 

15-by-30 

0.8526E+1 

0.4655E-5 

A2 

15-by-30 

0.1502E+1 

0.5494E-5 

A3 

15-by-30 

0 . 1069E+1 

0.0 

a4 

15-by-30 

0.2373E+1 

0.7918E-5 

A5(l,2,2) 

15-by- 30 

0.1108E+10 

0.0 

A5(l,3,4) 

15-by- 30 

0.9105E+14 

0.1304E-5 

Ag(l,l,2) 

15-by- 30 

0 . 5697E+I 

0.1931E-5 

AgUdd) 

15-by-30 

0. 4282E+3 

0.5759E-5 

A1 

15-by-45 

0.1554E+2 

0.4655E-5 

A2 

15_by-i)5 

0.2586E+1 

0.5494E-5 

A3 

15_by-H5 

0.1999E+1 

0.0 

a4 

15-by-45 

0 . 4002E+1 

0.7918E-5 

Aj(1,2,2) 

15-by-45 

0. 1686E+10 

0.0 

A5(l,3d) 

15-by-45 

0.1411E+15 

0.1304E-5 

Ag ( 1 > 1 » 2  ) 

15-by-45 

0.8284E+1 

0.1931E-5 

AgCl^,1*) 

15-by-45 

0.5404E+3 

0.5759E-5 

A1 

20-by-40 

0.1029E+4 

0.7707E-5 

A2 

20-by-40 

0. 3248E+3 

0.1008E-4 

A3 

20-by-40 

0 . 9273E+2 

0 . 2363E-6 
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TABLE  5.1  (Continued) 


MATRIX 

DIMENSION 

PREDICTED 
RELATIVE  ERROR 

ACTUAL 

RELATIVE  ERROR 

A4 

20-by-40 

0.4909E+3 

0.1090E-4 

A5(l,2,3) 

20-by-40 

0.1114E+15 

o 

• 

o 

A5(l,3,3) 

20-by-40 

0.1398E+21 

0.1232E+1 

a6(i,1,3) 

20-by-40 

0.1366E+4 

0 . 1120E-4 

Ag  ( 1 , 3 , 3 ) 

20-by-40 

0.7873E+5 

0.5103E-5 

> 
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5.2.2  Results  for  the  Pseudoinverse  Table  5.2  gives 
the  predicted  relative  error  and  the  actual  relative  error 
incurred  in  the  computation  of  the  Pseudoinverse  A+ , 
given  by  (5.19)  and  (5.25),  of  a  matrix  A  of  dimension 
m-by-2m  or  m-by-3m.  The  predicted  absolute  error, 
llAg-A+fl^,  is  given  by  (4.59),  and  the  predicted  relative 
error  is  defined  to  be 


The  actual  relative  error  can  be  computed  since  the  exact 
Pseudoinverse  is  known. 

The  first  column  of  Table  5.2  lists  the  nonsingular 
matrices  A.  of  Subsection  5.1.1  used  to  construct  the 

l 

matrices  A,  which  take  either  the  form  of  (5.14)  or 
(5.15).  The  three  parenthesized  numbers  after  A^  and 
Ag  are,  respectively,  the  values  of  a,  h  and  r  used 
to  construct  these  matrices.  The  second  column  of  the 
table  gives  the  dimension  of  A.  In  many  cases,  the 
condition  specified  by  (4.47)  does  not  hold,  and  no  pre¬ 
dicted  relative  error  can  be  computed  for  these  matrices. 
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TABLE  5.2 

ERROR  IN  THE  PSEUDOINVERSE 


MATRIX 

DIMENSION 

PREDICTED 
RELATIVE  ERROR 

ACTUAL 

RELATIVE  ERROR 

Ai 

3-by-6 

0.8974E-4 

0.2514E-5 

A2 

3-by-6 

0 . 1387E-1 

0.4023E-5 

A3 

3-by-6 

0 . 4948E-2 

0.1252E-5 

A4 

3-by-6 

0.1105 

0.3219E-5 

Ag(l, 1.5,1) 

3-by-6 

0.1595E-2 

0.5861E-5 

l - 1 

C\J 

*\ 

1 — 1 

l n 
< 

3-by-6 

0.3203E-3 

0 . 6606E-6 

Ag (1,2,1) 

3-by-6 

0 . 2692E-3 

0.9656E-6 

Ag(l,3,l) 

3-by-6 

0.1928E-3 

0.1460E-5 

A1 

3-by-9 

0.8974E-4 

0.2095E-5 

A2 

3-by-9  . 

0.1387E-1 

0 . 8404E-5 

A3 

3-by-9 

0.4948E-2 

0.8762E-5 

A4 

3-by-9 

0.1105 

0.1878E-5 

Ag  (1,1. 5,1) 

3-by-9 

0 . 1595E-2 

0. 3100E-5 

1 - 1 

r 

CM 

1 — 1 

LPv 

< 

3-by-9 

0.3203E-3 

0 . 4433E-6 

Ag (1,2,1) 

3-by-9 

0.2692E-3 

0.7544E-6 

Ag(l,3,l) 

3-by-9 

0.1928E-3 

0.9835E-6 

A1 

5-by-10 

0 . 5797E-2 

0.2490E-5 

A2 

5-by-10 

— 

0 . 3100E-4 

A3 

5-by-10 

0.4170 

0.9749E-5 

a4 

5-by-10 

— 

0.7093E-5 

•  , 
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TABLE  5.2  (Continued) 


MATRIX 

DIMENSION 

PREDICTED 
RELATIVE  ERROR 

ACTUAL 

RELATIVE  ERROR 

A5(l, 1.5,1) 

5-by-10 

0.4170E-2 

0.1171E-5 

(1,2,1) 

5-by-10 

0.9009E-3 

0.1008E-5 

Ag (1,2,1) 

5-by-10 

0.4285E-2 

0.3193E-5 

A6(i,3,l) 

5-by-10 

0. 3244E-2 

0.2076E-5 

A1 

5-by-15 

0.5797E-2 

0.2075E-5 

A2 

5-by-15 

— 

0.1886E-4 

A3 

5-by-15 

0.4171 

0 . 9521E-4 

A4 

5-by-15 

— 

0 . 1010E-4 

A5(l, 1.5,1) 

5-by-15 

0.4170E-2 

0.1794E-5 

1 - 1 

C\J 

•N 

1 - 1 

in 

< 

5-by-15 

0.9009E-3 

0.8430E-6 

Ag (1,2,1) 

5-by-15 

0. 4285E-2 

0.2638E-5 

Ag( 1,3,1) 

5-by-15 

0.3244E-2 

0.2103E-5 

A1 

8-by-l6 

— 

0.6243E-5 

A2 

8-by-l6 

— 

0.1755E-3 

A3 

8-by-l6 

— 

0 . 5575E-4 

A4 

8-by-l6 

— 

0.1529E-4 

A5(l, 1.5,1) 

8-by-l6 

0 . 1814E-1 

0 .  A313E-5 

(1,2,1) 

8-by-l6 

0.2918E-2 

0.1650E-5 

Ag (1,2,1) 

8-by-l6 

0.1125E+1 

0.5289E-5 

Ag(l,3,l) 

8-by-l6 

0.7700 

0.1771E-5 

A1 

8-by-24 

— 

0.4962E-5 

A2 

8-by-24 

— 

0.1630E-3 

(x«s,j;)aA 
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TABLE  5.2  (Continued) 


MATRIX 

DIMENSION 

PREDICTED 
RELATIVE  ERROR 

ACTUAL 

RELATIVE  ERROR 

A3 

8-by-24 

— 

0.6708E-3 

A4 

8-by-24 

— 

0.2341E-4 

A-Cl, 1.5,1) 

8-by-24 

0 . 1814E-1 

0.3664E-5 

Ag (1,2,1) 

8-by-24 

0.2918E-2 

0.2239E-5 

1 - 1 

C\J 

i — 1 

VO 

< 

8-by-24 

0.1125E+1 

0.1757E-5 

Ag(l,3,D 

8-by-24 

0.7700 

0 . 76-43E-5 

A1 

10-by-20 

— 

0.1233E-4 

A2 

10-by-20 

— 

0.2342E-3 

A3 

10-by-20 

— 

0.1702E-3 

a4 

10-by-20 

— 

0.4673E-4 

A5(l, 1.5,1) 

10-by-20 

0.5584E-1 

0.7403E-5 

Aj  (1,2,1) 

10-by-20 

0 . 5352E-2 

0.2238E-5 

Ag  (1,2,1) 

10-by-20 

— 

0.2979E-5 

Ag(l,3,l) 

10-by-20 

— 

0.1031E-4 

A1 

10-by-30 

— 

0.1078E-4 

A2 

10-by-30 

— 

o.is^e-s 

A3 

10-by-30 

— 

0.1699E-2 

a4 

10-by-30 

— — 

0.6107E-4 

A5(l, 1.5,1) 

10-by-30 

0.5584E-1 

0.731*9E-5 

Ag  (1,2,1) 

10-by-30 

0 . 5352E-2 

0 . 2367E-5 

Ag  (1,2,1) 

10-by-30 

— 

0.2339E-4 

Ag(l,3,l) 

10-by-30 

— 

0.1283E-4 

. 

110 


TABLE  5.2  (Continued) 


MATRIX 

DIMENSION 

PREDICTED 
RELATIVE  ERROR 

ACTUAL 

RELATIVE  ERROR 

Ai 

15-by-30 

— 

0.8180E-5 

A2 

15-by- 30 

— 

0.1050E-2 

A3 

15-by-30 

— 

0.1160E-2 

a4 

15-by-30 

— 

0.9946E-3 

a5(i, 1.5,1) 

15-by-30 

0.2512E+1 

0.5378E-5 

1 - 1 

C\J 

1 — 1 

LH 

< 

15-by-30 

0.1685E-1 

0.7805E-5 

Ag (1,2,1) 

15-by-30 

— 

0.2690E-5 

a6(i,3,d 

15-by-30 

— 

0.8783E-5 

A1 

15-by-45 

— 

0.6357E-5 

A2 

15-by-45 

— 

0.1738E-2 

A3 

15-by-45 

— 

0 . 1114E-1 

A4 

15-by-45 

— 

0 . 4640E-3 

A5(l, 1.5,1) 

15-by-45 

0.2512E+1 

0.4315E-5 

1 - 1 

C\J 

1 — 1 

l n 
< 

15-by-45 

0.1685E-1 

0.5256E-5 

Ag (1,2,1) 

15-by-45 

— 

0.1582E-4 

Ag(l,3,l) 

15-by-45 

— 

0.1282E-4 

A1 

20-by-40 

— 

0.1536E-4 

A2 

20-by-40 

— 

0 . 8582E-2 

A3 

20-by-40 

— 

0.3796E-2 

A4 

20-by-40 

— 

0 . 1100E-2 

A5(l, 1.5,1) 

20-by-40 

— 

0.2162E-4 

■ 
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TABLE  5.2 

( Continued ) 

MATRIX 

DIMENSION 

PREDICTED 
RELATIVE  ERROR 

ACTUAL 

RELATIVE  ERROR 

i — 1 

C\J 

1 — 1 

LP\ 

20-by-40 

0. 3908E-1 

0. 3918E-5 

i — 1 

C\ J 

i — 1 

VO 

20-by-40 

— 

0.2529E-4 

g ( 1 >  3 j 1) 

20-by-40 

— 

0.1741E-4 

CHAPTER  VI 


CONCLUSIONS  AND  SUGGESTIONS 
FOR  FUTURE  RESEARCH 


6 . 1  Conclusions 

The  following  conclusions  are  drawn  from  the  numerical 
results : 

1.  The  predicted  round-off  error  unfortunately  is 
substantially  larger  than  the  actual  round-off  error  in 
most  cases.  This  is  due  to  the  fact  that  in  the  error 
analysis,  we  allow  for  the  maximum  amount  of  round-off 
error  in  each  calculation.  In  practice,  the  round-off 
error  incurred  in  a  computation  is  seldom  the  maximum 
possible.  It  is  important,  however,  to  realize  that 
matrices  can  be  constructed  such  that  the  actual  total 
round-off  error  approaches  the  predicted  value. 

2.  The  comparatively  large  values  obtained  for  the 

predicted  relative  error  in  computing  the  Reflexive 

Generalized  Inverse  of  a  matrix  constructed  from  Ac  are 

b 

due  to  large  values  of  the  bound  on  flF.J^.  Consequently, 
the  computation  of  (A  ) ’  should  be  done  in  multiple 
precision  arithmetic  to  lower  the  bound  on  ||  F ^ ||  ^ ,  and 
thus  improve  the  bound  for  the  Reflexive  Generalized 


Inverse . 


« 
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3.  The  actual  relative  error  is,  on  the  average, 
smaller  for  the  computed  Reflexive  Generalized  Inverse, 
indicating  that  its  computation  is  less  susceptible  to 
the  accumulation  of  round-off  errors  than  is  the  computa¬ 
tion  of  the  Pseudoinverse.  This  remark,  however,  applies 
only  to  the  two  algorithms  that  we  have  programmed,  and  is 
based  solely  on  the  set  of  data  tested. 

4.  The  complete  listing  of  the  numerical  results 
shows  that  the  bound  given  by  (4.59)  is  approximately 
equal  to  the  product  of  11?*!^  and  II  ^4 II  ^  •  Hence,  these 
are  the  only  two  values  which  need  be  computed  to  obtain 

an  approximate  bound  on  the  round-off  error  in  the  computed 
Pseudoinverse . 

6 . 2  Suggestions  for  Future  Research 

1.  In  some  applications  of  generalized  inverses,  it 
is  irrelevant  which  of  the  four  generalized  inverses 
defined  in  Chapter  I  is  computed.  Hence,  error  analyses 
of  other  algorithms  could  suggest  which  of  the  four 
generalized  inverses,  and,  in  particular,  which  algorithms, 
are  the  least  susceptible  to  the  accumulation  of  round-off 
error . 

2.  The  error  bound  for  the  matrix  B#,  ||  II 
Winner’s  Algorithm  would  be  significantly  smaller  for 
large  matrices  if  the  “-norm  replaced  the  1-norm  in  (4.27). 


« 


. 
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This  would  require  use  of  the  1-norm,  instead  of  the 
°°-norm,  in  bounding  the  error  for  the  computed  Hermite 
normal  form.  The  error  bound  for  the  computed  Pseudo¬ 
inverse  may  then  be  significantly  smaller,  because,  from 
observation  of  the  numerical  results,  ||  P||  ,  which  is  a 

factor  in  the  bound  for  IlF^II*,,  is  approximately  equal 
to  II  P*||  ooll  A|l  a,!!  oo  •  As  stated  in  the  fourth  conclusion 
of  Section  6.1,  HFjJ^  is  one  of  the  two  major  factors 
of  the  bound  on  the  computed  Pseudoinverse. 

3.  A  major  deficiency  in  the  bound  given  by  (4.59) 
for  the  round-off  error  incurred  in  computing  the  Pseudo¬ 
inverse  is  the  requirement  that  condition  (4.47)  hold. 
This  condition  might  be  relaxed  if  a  simple  relationship 
could  be  found  between  the  exact  and  computed  values  of 
||  (PAB+F)-1||  .  Perhaps,  for  example,  it  could  be  shown 

that 


||  (PAB+F)"1!^  <  ||  (PAB+F)"1||oo, 

provided  the  condition  number  of  (PAB+P)  is  greater  than 
or  equal  to  one. 

4.  The  predicted  relative  error  for  the  computed 
Reflexive  Generalized  Inverse  varies  greatly  with  different 
test  matrices  of  the  same  dimension,  although  the  actual 
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relative  error  is  quite  stable.  Perhaps  large  values  of 
the  predicted  relative  error  could  be  related  to  a  con¬ 
dition  number  or  a  norm  of  the  matrix.  This  would  give 
a  criterion  for  deciding  whether  or  not  the  predicted 
relative  error  is  a  reasonable  approximation  to  the 
actual  relative  error. 
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APPENDIX 


LISTINGS  OF  FORTRAN  IV  SUBPROGRAMS 

The  following  pages  contain  listings  of  the  IBM  System 

360  FORTRAN  IV  subprograms  used  to  compute  a  Reflexive 

Generalized  Inverse  and  the  Pseudoinverse  of  a  matrix  A. 

The  SUBROUTINE  AGI  computes  a  Reflexive  Generalized  Inverse, 

and  the  SUBROUTINE  PSEUDO  computes  the  Pseudoinverse.  The 

SUBROUTINE  NORM  computes  the  “-norm  of  an  arbitrary  matrix. 

The  FUNCTION  subprogram  AMAX  computes  the  element  of  maximum 

modulus  in  a  certain  submatrix  of  an  arbitrary  matrix. 

The  unit  round-off  error  for  the  IBM  System  360/67 
-21 

computer  is  2  .  Although  the  word  length  is  32  bits,  one 

bit  is  used  for  the  sign,  seven  bits  are  used  for  the 
exponent,  and  the  three  leading  bits  of  the  fractional  part 
may  be  zero,  as  normalization  is  performed  in  hexadecimal. 

The  matrices  of  round-off  errors  are  referenced  in 
comment  statements  by  the  names  given  them  in  Chapters  III 
and  IV. 
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SUBROUTINE  AG  I < M , N , ZERO , A, Q, GNORM , R EL ERR ) 

0 1  ME NS  ION  A ( 50,50) ,Ul 50,50) ,P( 50,50) ,PD( 50,50 ) , 

1  Q(5C,50),NR<50),NC<50) ,DEL(50) 

DOUBLE  PRECISION  E , EP S I , T EMP , F2NORM , E 1 , E2 , DP , DM ( 5 1 ) 
C 

C  THIS  SUBROUTINE  CALCULATES  A  REFLEXIVE  GENERALIZED 
C  INVERSE  OF  A  MATRIX  A  OF  DIMENSION  M-BY-N  AND 
C  CALCULATES  AN  UPPER  BOUND  FOR  THE  ROUND-OFF  ERROR 
C  INCURRED  IN  THE  CALCULATION*  THE  REFLEXIVE  GENERALIZED 
C  INVERSE  IS  RETURNED  VIA  THE  ARRAY  U. 

C  GNORM  IS  THE  VALUE  OF  THE  INFINITY  NCRM  OF*  THE  EXACT 
C  REFLEXIVE  GENERALIZED  INVERSE  OF  A.  IT  IS  USED  TO 
C  COMPUTE  THE  PRELICTED  RELATIVE  ERROR,  RELERR,  IN  THE 
C  COMPUTED  REFLEXIVE  GENERALIZED  INVERSE. 

C  EPSI  IS  THE  UNIT  ROUNO-OFF  ERROR  FOR  THE  IBM  360/b7. 

C  ZERO  IS  THE  ZERO  CRITERION  USED  IN  THE  FORWARD 
C  ELIMINATION. 

C 

EPSI=2.D0**(-21) 

FORM  A  UNIT  MATRIX  IN  U. 

DATA  U/25C0*C/ 

DO  10  1=1,50 
10  U(I,I)=1 

NRANK  IS  THE  CALCULATED  RANK  OF  A.  G  IS  THE  VALUE  OF 
THE  ELEMENT  OF  MAXIMUM  MODULUS  IN  THE  REDUCTION  OF  A. 

MM 1=M— 1 
MP 1=M+1 
MP2=M+2 
NRANK=M 
NFLAG=0 
G=0 

THE  BEGINNING  OF  THE  FORWARD  ELIMINATION. 

OC  45  1=1,  M 

G=AMAX1(G,AMAX( I , M , N , A ) ) 

IP1=I+1 

THE  COMPLETE  PIVOTING.  THE  VECTORS  NR  ANO  NC  NOTE 
THE  NECESSARY  ROW  AND  COLUMN  INTERCHANGES  RES PEC  I T V EL Y . 


DMAX=0 
NR ( I )  =  I 
NC ( I ) = I 


. 
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DO  20  J=I,M 
DO  20  K=  I ,  N 

IFIDMAX.GE.ABSI A( J,K) I)  GOTO  20 
DMAX=ABS(AIJ,K) ) 

NR  I  I ) = J 
NCII)=K 
20  CONTINUE 

IFIDMAX.LT. ZERO)  NRANK= 1-1 
C 

C  IF,  AT  ANY  STAGE  OF  THE  ELIMINATION,  THE  REMAINDER  OF 
C  THE  MATRIX  IS  LESS  THAN  THE  ZERO  CRITERION,  THEN  NO 
C  FURTHER  ELIMINATION  IS  NECESSARY,  AND  THE  RANK  UF  A  IS 
C  FOUND  TO  BE  LESS  THAN  M. 

C 

IFINRANK.NE.M)  GCTO  50 

INTERCHANGE  COLUMNS  IF  NECESSARY. 

1FINCII  J.EQ.I)  GOTO  30 
DO  25  K= 1 , M 
DMAX= A  I K  ,  I  ) 

AIKfl )=A(K,NC(I)) 

25  AIK, NCI  I ) )  =  DMAX 

30  IFII.EQ.M)  GOTO  50 

INTERCHANGE  ROWS  IF  NECESSARY. 

IFINRID.EO.I)  GOTO  36 
DC  35  J=I ,N 
LMAX=A I  I , J) 

A 1 1 ,  J  )— A  I  NR  I  I J  ,  J  ) 

35  AINRI I ) , J ) -UMAX 

DETERMINE  WHETHER  OR  NOT  THE  REMAINDER  OF  THE  MATRIX 
IS  ALL  ZEROS. 

36  CMAX=0 
CO  40  J-  I  P 1 ,  M 
DC  40  K= I , N 

IFIDMAX.GE.ABSI A( J , K ) ) )  GOTO  40 
DMAX= ABS I  A I J  ,  K ) ) 

40  CONTINUE 

iFIDMAX.LT  .ZERO)  NR  ANK-=  I 
IFINRANK.NE.M)  GOTO  50 

THE  REDUCTION  OF  ROWS  <1*1,  ...»  M)  OF  A. 


. 
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00  45  I  P  1  f  M 
CP=-A  (J«I)/MI«I) 

A(  J,  I  )=0P 
DO  45  K=IP1,N 

45  A< JtK)=A(J,K)+DP*A( I ,K) 
INITIALIZE  THE  MATRIX  P. 


50 

lF(NRANK.NE.Pl) 
NRANK  =  MM  1 
NFLAG= 1 

GOTO  51 

51 

DG  65  J= 1 , K 

DO  65  K=  1 1  M 

65 

P( J,K)=U( J,K) 
NRP1=NRANK*  1 

DC  68  I=NRP1,P 

68 

P( I «  NRANK )  =  A  (  I 

, NRANK) 

CALCULATION  OF  P. 

00  55  1= 1 , NRANK 
ID=NRP1-I 


INTERCHANGE  COLUMNS  IF  NECESSARY  • 

I F  ( N  K I  ID)  .EQ. ID)  GOTO  71 
00  70  J=1,M 
DPAX=P( J,  ID) 

P( Jt ID)=P( JtNR(ID) ) 

7C  P  (  J  »  NR ( ID) ) -DP  AX 

71  IF( I .EQ. NRANK)  GOTO  95 

FORM  THE  MATRIX  P(IC-l)  IN  PU. 

DO  75  J=1»M 
DO  75  K=  1 ,  P 
75  P0( JtK)*U( JiK) 

DO  80  J=  1 0  f  M 

80  PD(J,  ID-1 )  =  A ( J  *  1 0- 1 ) 

MULTIPLY  P  BY  PC. 

00  90  J=  1 ,  P 
DO  85  K*ltP 
OP ( K ) *0 
DO  85  L-  1 » M 

85  OMiK)  =  DM(K)+DBLE(PUtL)  )*PD(L,K) 
DC  90  K*ltP 
90  P(JtK)-DMK) 

95  CONTINUE 


n  n  r  n  o  o  o  noo  oooooo 


THE  INFINITY  NORM  OF  P  IS  STORED  IN  PNORM. 

CALL  NORM(PNORM,M,M,P) 

CALCULATION  CF  C. 

NRANK=NR ANK+NFL AG 
NRP 1=NRANK+ 1 
NRP2=NR ANK+2 
DO  IOC  1=1, N 
DG  100  J= 1 , N 
100  0<I,J)  =  UU,J) 

OC  106  1  =  1,  NR ANK 
IP1=I+1 

DC  106  J= I P  1 ,  NRP 1 
L=  J- 1 
DP  =  C 

DC  105  K= I , L 

105  DP=DP-UbLE(A(K,J) )/A(K,K)*Q(I,K) 

106  Q ( I , J )=DP 

I F ( N  •  LE • NRP 1 )  GOTO  112 
DO  111  1  =  1 , NR ANK 
DO  ill  J=NRP2,N 
DP  =  0 

DO  11C  K=I , NR ANK 

110  DP  =  DP-DBLt(MK,J)  )  /  A  (  K  ,  K )  *0  (  I  ,  K ) 

111  Q ( I , J ) =DP 

INTERCHANGE  THE  RCwS  OF  g  IF  NECESSARY. 

112  DO  1 2  C  1  =  1 ,  NR ANK 
I  D=NRP 1-1 

If (NC(ID). EC. IG)  GOTO  120 
DC  115  J= 1 , N 
GMAX=C ( NC ( IC) , J) 
g(NC(  ID) • J)=ti(  ID, J) 

115  Q  ( I  G , J )  =  CP AX 
120  CONTINUE 

ThE  INFINITY  NORP  CF  0  IS  STORED  IN  QNORM. 

CALL  NCRMUNCRP,N,N,0) 

CALCULATION  CF  THE  PATRIX  DELTA-.  ITS  “DIAGONAL 
ELEMENTS  ARE  STCREC  IN  THE  VECTOR  DEL. 

DO  125  1  =  1,  NR ANK 
125  DEL (  I  )= 1/ A (  I , I ) 

IF(NRANK.EC.P)  GOTO  131 


. 
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DC  13C  I=NRP1,M 

130  DfcL (  I  )  =  0 

THE  INFINITY  NORM  OF  DELTA-  IS  STORED  IN  DELNRM • 

131  DtLNRM=0 

DO  135  I  =  1 »  NK ANK 

135  DELNRM=AMAX1  (DELNRM,  ABS  (  DEL  (  1)11 
ThE  INFINITY  NORN  GF  THE  MATRIX  F2  IS  STOR-ED  IN  F2NURM. 
F  2NLRM  =  0 

DO  140  I=1,NRANK 
fc= ( I-1)*2.C1*G*EPSI 

T EMP=DABS (  (DELI  I ) *EPSI-E )/ < CEL ( I  )*(DEL(  I )+E) ) ) 

140  F2NORM=DMAXl(F2NORM,TEMP) 

CALCULATION  OF  THE  PREDICTED  ABSOLUTE  ERROR  IN  THE 
COMPUTED  REFLEXIVE  GENERALIZED  INVERSE. 

E1=NRANK*( N-NRANK+1 ) *2** < NRANK-2 ) *EPSI 
£2=DELNRM+F2NCRM 

ERROR* <El*fc2  +  CNCRM*F2N0RM)*PNGRM  +  ( NR ANN- 1  )  * 

1  2**NRANK*EPSI*E2*(QN0RM*E1 )+2*EPSI*0N0KM* 

2  D£LNRM*PNORM 
RELERR=ERRCR/GNGRM 


MULTIPLICATION  CF  Q  eY  DELTA-. 

DC  145  1*1,  M 
DO  145  J=1,N 

145  U(J,IJ  =  Q(Jf  n*DEL(I) 

MULTIPLICATION  CF  ( C*DELT A- )  BY  P.  THE  REFLEXIVE 
GENERALIZED  INVERSE  CF  A  IS  STORED  IN  THE  FIRST  M 
CGLUMNS  CF  Q. 

DO  155  I  =  1 1  N 
DC  150  J  =  1 ,  M 
D  M  (  J  )  =  0 
DC  15C  K*1,M 

15C  DM(J)-DM(J)+DbLE(C(I»K) ) *P ( K , J ) 

DC  155  J=  1 »  M 
155  (.(I,  J)*DM(  J) 

RETURN 

END 


■ 
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SUBROUTINE  PS EUOOi M , N, Z ERO , A, G, GNORM , RELERR ) 
DIMENSION  A ( 50 »50) , I DENT ( 50) , AA(  50, 50) , P ( 50, 50 ) , 

1  6(50,50) 

DOUBLE  PRECISION  DP,DM,EPSI 
C 

C  THIS  SUBROUTINE  CALCULATES  THE  PSEUDO  I NVERSE  OF  A 
C  MATRIX  A  OF  DIMENSION  M-BY-N  AND  CALCULATES  AN  UPPER 
C  BOUND  FOR  THE  ROUND-OFF  ERROR  INCURRED  IN  THE 
C  COMPUTATIONS*  THE  PSEUDO  I NVERSE  IS  RETURNED  VIA  THE 
C  ARRAY  G. 

C  ZERO  IS  THE  ZERO  CRITERION  USED  IN  THE  FORWARD 
C  ELIMINATION.  GNORM  IS  THE  VALUE  OF  THE  INFINITY  NORM 
C  OF  THE  EXACT  PS EUDOI NVERSE  OF  A.  IT  IS  USED  TO  COMPUTE 
C  THE  PREDICTED  RELATIVE  ERROR,  RELERR,  IN  THE  COMPUTED 
C  PSEUDOINVERSE. 

C  EPSI  IS  THE  UNIT  ROUND-OFF  ERROR  FOR  THE  IBM  360/67. 

C 

EPSI=2.D0**(-21) 

THE  MATRIX  A  IS  STORED  IN  AA  FOR  LATER  USE. 

THE  INFINITY  NORM  OF  A  IS  STORED  IN  ANORM. 

DO  10  1=1, M 
DO  10  J=  1 ,  N 
10  AA( I , J)=A( I , JJ 

CALL  NGRM( ANORM,M,N, A) 

AT  EACH  STAGE  OF  THE  REDUCTION,  NC  IS  THE  NUMBER  OF  THE 
COLUMN  OF  A  IN  WHICH  THE  PIVOTAL  ELEMENT  IS  LOCATED. 
THUS  COLUMN  NC  OF  A  WILL  BE  REDUCED  TO  A  COLUMN  VECTOR 
OF  THE  UNIT  MATRIX  AT  EACH  STAGE. 


NC=0 


THE  BEGINNING  OF  THE  CALCULATION  OF  THE  HERMI TE  NORMAL 
FORM  OF  A.  F2N0RM  DENOTES  THE  INFINITY  NORM  OF  THE 
MATRIX  OF  ROUND-OFF  ERRORS  INCURRED  IN  COMPUTING  THE 
HERMITE  NORMAL  FORM  OF  A. 

F2NQRM=1 
00  70  1=1, M 
NC=NC+ 1 

IF(NC.GT.N)  GOTO  75 

NCP 1=NC+ 1 

IP1=1+1 

USE  PARTIAL  PIVOTING  TO  FIND  THE  PIVOTAL  ELEMENT  FOR 
THE  I-TH  ROW. 


CMAX=ABS(A(i,NC)) 


•  ) 
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L=  I 

IFll.fcg.M  GOTO  21 
DG  20  J=  I  P  1 ,  M 

IF(DMAX.GE.AbS(A(J,I) ))  GCTU  20 
L  =  J 

UMAX= A8S ( A  (  J  ,  I ) ) 

20  CONTINUE 

21  IF(DMAX.GT.ZERG)  GOTO  40 
IFUC.EG.NJ  GOTO  75 

IF  ALL  GF  THE  POSSIBLE  PIVOTAL  ELEMENTS  IN  COLUMN  NC 
ARE  ZERO,  SEARCH  FOR  A  PIVOTAL  ELEMENT  IN  COLUMNS 

(NC+1,  ...  ,  N). 

DO  30  J  =  NCP  1 ,  N 
NC=NC+ 1 

USE  PARTIAL  PIVOTING  TO  LUCATL  THE  LARGEST  POSSIBLE 
PIVOTAL  ELEMENT  IN  COLUMN  J  • 

DM AX= A8S ( A (  I,J) ) 

L=  I 

I F ( I .EQ.M)  GOTO  26 
DO  25  K=  I  P  1 ,  M 

I F ( DMAX . GE • A6S i A(K, J) ) )  GOTO  25 
L=K 

DMAX= A8S ( A  (  K  ,  J  )  I 

25  CONTINUE 

26  NCP 1=NC+ 1 
IF(DMAX.GT.ZERC)  GOTO  40 

30  CONTINUE 

AS  ALL  CF  THE  POSSIBLE  PIVOTS  ARE  ZERO,  THE  HtRMITE 
NORMAL  FORM  CF  A  HAS  BEEN  FOUND • 

GOTO  75 

40  IF(L.EQ.I)  GOTO  46 

INTERCHANGE  ROUS  IF  NECESSARY. 

DO  45  J=NC,N 
CMAX= A ( I  ,  J) 

AUt  J)*A(Lf  J) 

45  A ( L , J )-UM AX 

46  IF(NC.EU.N)  GOTO  61 

THE  REDUCTION  GF  COLUMNS  (NC+1,  ...  ,  Nl  OF  A  USING 
THE  FORWARD  ELIMINATION  OF  GAUSSIAN  ELIMINATION  WITH 
ACCUMULATION  OF  INNER  PRODUCTS. 


. 
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CM=0 

DC  55  J=1,M 
IF(J.EQ.I)  GOTO  55 
DP=— A(J,NC)/A( ItNCI 
DM=DM+ABS(DP) 

DO  50  L—NCP  1 »  N 
5C  A(  J,L)=A(  J,L  )+DP*A<  I  ,L) 

55  CONTINUE 

DM=0M+1/ABS(A(I,NC>) 

F2NCRM=F2N0RM*DMAX1( 1.00, DM) 

DIVIDE  THE  I-TH  ROW  OF  A  BY  THE  PIVOTAL  ELEMENT. 

DO  60  J=NCP 1 « N 

60  A<I,J)=A(I, J)/A<I,NC) 

SET  COLUMN  NC  OF  A  EQUAL  TO  THE  I-TH  COLUMN  OF  THE  UNIT 

MATRIX. 

61  DO  65  J*1#M 

65  A( J  t  NC 1-0 

All  ,NC)®1 

ThE  VECTOR  IOENT  INDICATES  WHICH  COLUMNS  OF  THE  HERMITE 
NORMAL  FORM  OF  A  ARE  EQUAL  TO  COLUMNS  OF  THE  UNIT 
MATRIX.  NRANK  IS  THE  RANK  OF  A. 

IDENT ( I j=NC 
70  NRANK= I 

75  F2NGRM=NRANK*EPSI*F2N0RM 

CALCULATION  OF  THE  MATRIX  P-TRANSPOSE. 

DO  80  1  =  1 1 NRANK 
J=IDENT ( I ) 

DO  80  K=  1 ,  M 
80  P ( 1 1 K )  =  AA ( K  » J ) 

THE  INFINITY  NORM  OF  P-TRANSPOSE  IS  STORED  IN  PNORM. 

CALL  NORMIPNCRM, NRANK, M,P) 

MULTIPLY  P-TRANSPOSE  BY  A  USING  ACCUMULATION  OF 
INNER  PRODUCTS. 

DO  86  1=1, NRANK 
DO  86  J= 1 , N 
D  P=0 

DO  85  K= 1 , M 

0P=DP+DBLE4P< I ,KJ ) *AA(  K, J) 


85 


o  o  o  n  o  oooo  non  oooooo  o  o  o  o  ooooo 


129 


86  G  ( 1 ,  J )  =  DP 

THE  INFINITY  NORM  OF  P-TRANSPOSE  *  A  IS  STURED  IN 
PANORM.  F1N0RM  IS  THE  INFINITY  NORM  OF  THE  MATRIX  OF 
ROUND-OFF  ERRORS  INCURRED  IN  COMPUTING  P-TRANSPOSt  *  A. 

CALL  NCRMI PANORM , NR ANK , N , G ) 

F1NCRM=EPSI*PNQRM*AN0RM 

MULTIPLY  (P-TRANSPOSE  *  A)  BY  B-TRANSPOSE  USING 
ACCUMUL AT I  ON  OF  INNER  PRODUCTS. 

DO  91  1= 1  *  NR ANK 
DO  91  J=  1 ,  NRANK 
0  P  =  Q 

DO  90  K=  1 ,  N 

90  DP=DP  *DBLc ( G ( I , K ) )*A( J,K) 

91  AA (  I  ,  J )  =DP 

INVERT  THE  MATRIX  (P-TRANSPCSE  *  A  *  B-TR AN SPCSt ) . 

CALL  THE  INVERTED  MATRIX  Y. 

GG  DENOTES  THE  ELEMENT  OF  MAXIMUM  MUOULUS  IN  THE 
COMPUTATION  OF  Y. 

CALL  A INV  ( AAyGfNRANKfGG) 

THE  INFINITY  NCRM  OF  Y  IS  STORED  IN  PABNRM . 

CALL  NCRM(PABNRM, NRANK, NRANK, G) 

MULTIPLY  B-TRANSPOSE  BY  Y  USING  ACCUMULATION  OF 
INNER  PRODUCTS. 

DO  96  1=1, N 
DO  96  J= 1 , NRANK 
DP=0 

DO  95  K=l, NRANK 

95  DP=DP+DBLE(A(K,  I) )*G(K, J) 

96  AA ( 1 , J )=DP 

MULTIPLY  (B-TRANSPOSE  *  Y)  BY  P-TRANSPOSE  USING 
ACCUMULATION  OF  INNER  PRODUCTS.  THE  RESULT,  THE 
PSEUDOINVERSE  OF  A,  IS  STORED  IN  G. 

DO  99  1=1, N 
DO  99  J= 1 , M 
DP=0 

DO  98  K=l, NRANK 

98  DP=DP+DBLEI AA ( I , K) )*P(K,J) 

99  G( I , J ) =DP 


' 


. 


o  o  o 
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THE  MATRIX  B-TRANSPOSE  IS  STORED  IN  P. 

DO  100  1  =  1*  NRANK 
DC  100  J= 1  *  N 

ioo  p(j»n=A(i»j) 
c 

C  THE  INFINITY  NORMS  OF  THE  MATRICES  B-TRANSPOSE,  F3,  F, 

C  INVERSE  OF  (P-TRANSPOSE  *  A  *  B-TRANSPOSE  +F ) ,  AND  F4 
C  ARE  STORED  IN  BNORM ,  F3N0RM ,  FNURM ,  PABFNM*  AND  F4N0RM 
C  RESPECTIVELY.  C  AND  0  ARE  CONSTANTS.  THE  PREDICTED 
C  VALUE  OF  THE  ABSOLUTE  ROUND-OFF  ERROR  IS  STORED  IN 
C  ERROR.  THE  PREDICTED  RELATIVE  ERROR  IS  STORED  IN 
C  RELERR. 

C 

CALL  NORM  CBNORM,N*NRANK*P) 

F3NGRM=EPSI*PAN0RM*8N0RM 

FN0RM=D6LE ( PNGRM) *ANGRM*F2N0RM+DBLE ( F 1N0RM ) *BNORM+ 

1  DBLE (F3N0RM ) 

C=GG*EPS I*NRANK*NRANK* ( 2 .005+NRANK ) 

TEST=1.D0-DBLE( FNORM ) *PABNRM 
D=TE  ST**2-4 .00*DBLE ( C ) *PAB  NRM 
PABFNM=.5CQ*(DBLE( TEST)-SQRT(D) )/C 
F4N0RM=DBLE ( P ABFNM ) * ( C*PAB FNM+FNORM*PABNRM ) 
ERROR=DBLEl PNORM)* (DBLE (BNORM ) *F4N0RM+DBL E < P ABNRM ) * 
1  ( F2N0RM+2*EPS I*BNQRM )  ) 

RELERR=ERROR/GNQRM 

RETURN 

END 


- 
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SUbRGUT I  NE  NGRM(X,MI,N,A) 

DIMENSION  A ( 50 , 50 ) 
oCUbLE  PRECISION  DP 

THIS  SUbRCJUT  i  NE  COMPUTES  THE  1NP1NITY  NORM  UP  1  HE 
MATRIX  A  ANL  STORES  IT  IN  X.  INNER  PRODUCTS  ARE 
ACCOMULATtU. 


X  =  0 

CG  10  1=1, M 
CP  =  0 

DO  5  J=1,N 

5  DP=LP+AbS l A l I  ,  J)  ) 
IP(LP-X)  1C, 1C, 6 

6  X  =  0  P 

10  CONTINUE 
RETURN 
END 
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FUNCTION  AMAX( I ,M,N, A) 

DIMENSION  A ( 50  »  50 ) 

THIS  FUNCTION  SUBPROGRAM  COMPUTES  THE  ELEMENT  OF 
MAXIMUM  MODULUS  IN  ROWS  I,...,M  AND  COLUMNS  Iy...,N  OF 
THE  RECTANCULAR  MATRIX  A  OF  DIMENSION  M-BY-N. 

AMAX=0 
DO  10  J=  1 1 M 
DO  10  K= I y  N 

10  AMAX=AMAX1(AMAX,ABS( A( J,K) ) ) 

RETURN 

END 


